How to use retime to get the plot of mean and stdev in continous time series
Show older comments
Hello there,
I have a datasets containing y-values in multi-years time series. My intention is getting the plot of monthly mean and its shaded error bars. The datasets contain long span of months which crossing the following year, i.e. in my datasets (attached) there from February, 2020 until December 2022.
But, at first, I want to get the data cleaned by selecting only y-value with a certain range. Here is my first sight code:
T = readtable('datamine');
Tr = table2timetable(T, 'RowTimes','time');
Tx = table(Tr.time,Tr.value_y,'VariableNames',{'time_a','val'});
idx = Tx.val>= 12 & Tx.val< 13.4; % selecting only y-value with a certain range
Tx_mean = retime(Tx(idx,:),'monthly','mean'); % get monthly averaged values ==> Failed
Tx_std = retime(Tx(idx,:),'monthly',@std); % get monthly stdev values ==> Failed
I tried to use directly retime function but it failed. Anyone know to get continous monthly mean and its shaded error bars properly?
More or less, the plot should be like this. The x-axis is continous months, crossing the following years:
Thanks!
Accepted Answer
The correct way to express the deviation of cata around the mean is to use the standard error of the mean, given by:
where σ is the standard deviation and N are the number of data used to calculate it.
T1 = readtable('datamine.xlsx')
TT1 = table2timetable(T1)
SEM = @(x) std(x)/sqrt(numel(x));
TT1momean = retime(TT1, 'monthly', 'mean')
TT1mosem = retime(TT1, 'monthly', SEM)
TT1monum = retime(TT1, 'monthly', 'count');
[Nmin,Nmax] = bounds(TT1monum{:,1})
% TT1Time = TT1mosem.time;
% TT1SEM = TT1mosem{:,1};
figure
plot(TT1momean.time, TT1momean{:,1}, '-k')
hold on
patch([TT1mosem.time; flip(TT1mosem.time)], [TT1momean{:,1}-TT1mosem{:,1}*1.96; flip(TT1momean{:,1}+TT1mosem{:,1}*1.96)], 'r', 'FaceAlpha',0.25, 'EdgeColor','r')
hold off
grid
xlabel('Time')
ylabel('Value')
title('Mean ±95% CI')
The SEM values are quite small when compared to the mean values (on the order of 
) so they are barely visible. I did a separate accumulation for the number of values in each month, and since they were all above 330, using the 95% confidence intervals from the normal distribution is a safe estimate. It is not necessary to use the t-distribution, since it will closely approximate the normal distribution with this many degrees-of-freedom.
) so they are barely visible. I did a separate accumulation for the number of values in each month, and since they were all above 330, using the 95% confidence intervals from the normal distribution is a safe estimate. It is not necessary to use the t-distribution, since it will closely approximate the normal distribution with this many degrees-of-freedom. .
12 Comments
Adi Purwandana
on 19 Oct 2024
Edited: Adi Purwandana
on 19 Oct 2024
Adi Purwandana
on 19 Oct 2024
As always, my pleasure!
My apologies — I forgot about the restriction.
I included it and the boxchart here —
T1 = readtable('datamine.xlsx')
idx = (T1.value_y >= 12) & (T1.value_y < 13.4); % selecting only y-value with a certain range
T1 = T1(idx,:)
TT1 = table2timetable(T1);
SEM = @(x) std(x)/sqrt(numel(x));
TT1momean = retime(TT1, 'monthly', 'mean');
TT1mosem = retime(TT1, 'monthly', SEM);
TT1monum = retime(TT1, 'monthly', 'count');
[Nmin,Nmax] = bounds(TT1monum{:,1})
% TT1Time = TT1mosem.time;
% TT1SEM = TT1mosem{:,1};
figure
plot(TT1momean.time, TT1momean{:,1}, '-k')
hold on
patch([TT1mosem.time; flip(TT1mosem.time)], [TT1momean{:,1}-TT1mosem{:,1}*1.96; flip(TT1momean{:,1}+TT1mosem{:,1}*1.96)], 'r', 'FaceAlpha',0.375, 'EdgeColor','r', 'EdgeAlpha',0.5)
hold off
grid
xlabel('Time')
ylabel('Value')
title('Mean ±95% CI')
[GV, ID1, ID2] = findgroups(year(T1.time), month(T1.time)); % Grouping Variable For ‘boxchart’
% [b1,b2] = bounds(ID1)
xtl = compose("%s",datetime(ID1,ID2,ones(size(ID1)), Format="MMM yyyy")); % Cell Array Of Time Values
% disp(xtl)
figure
hb = boxchart(GV, T1.value_y);
hb.JitterOutliers = "on";
hb.MarkerStyle = ".";
hb.Notch = "on";
hb.BoxMedianLineColor = "r";
Ax = gca;
xt = Ax.XTick;
Ax.XTick = [1 xt];
Ax.XTickLabel = xtl(Ax.XTick);
Ax.XTickLabelRotation = 45;
xlabel("Time")
ylabel("Value")
title("T1: Boxchart")
axis("padded")
% get(hb)
The boxchart took a bit of time because I do not use it frequently. It uses the original ‘T1’ with the ‘idx’ constraint, not any of the derived timetable arrays.
I defined the XTick locations as:
Ax.XTick = [1 xt];
to include thee first element. It needs to be integers from 1 to 35, within that constraint you can define it however you want. You can change the XTickLabel format by changing the Format string in the ‘xtl’ assignment. (That is the easiest way to do it in this instance.)
.
Adi Purwandana
on 19 Oct 2024
Great! But, I don't know it didn't work on my 2022a:
- the Xtick label shown as integer instead of MMM yyy
- the BoxMedianLineColor red doesn't show
Are the date time axis and BoxMedianLineColor working only on 2024b, as what you did?
clear all;
clc
T1 = readtable('datamine.xlsx')
idx = (T1.value_y >= 12) & (T1.value_y < 13.4); % selecting only y-value with a certain range
T1 = T1(idx,:)
TT1 = table2timetable(T1);
SEM = @(x) std(x)/sqrt(numel(x));
TT1momean = retime(TT1, 'monthly', 'mean');
TT1mosem = retime(TT1, 'monthly', SEM);
TT1monum = retime(TT1, 'monthly', 'count');
[Nmin,Nmax] = bounds(TT1monum{:,1})
% TT1Time = TT1mosem.time;
% TT1SEM = TT1mosem{:,1};
figure
plot(TT1momean.time, TT1momean{:,1}, '-k')
hold on
patch([TT1mosem.time; flip(TT1mosem.time)], [TT1momean{:,1}-TT1mosem{:,1}*1.96; flip(TT1momean{:,1}+TT1mosem{:,1}*1.96)], 'r', 'FaceAlpha',0.375, 'EdgeColor','r', 'EdgeAlpha',0.5)
hold off
grid
xlabel('Time')
ylabel('Value')
title('Mean ±95% CI')
[GV, ID1, ID2] = findgroups(year(T1.time), month(T1.time)); % Grouping Variable For ‘boxchart’
% [b1,b2] = bounds(ID1)
xtl = compose("%s",datetime(ID1,ID2,ones(size(ID1)), Format="MMM yyyy")); % Cell Array Of Time Values
% disp(xtl)
figure
hb = boxchart(GV, T1.value_y);
hb.JitterOutliers = "on";
hb.MarkerStyle = ".";
hb.Notch = "on";
hb.BoxMedianLineColor = "r";
Ax = gca;
xt = Ax.XTick;
Ax.XTick = [1 xt];
Ax.XTickLabel = xtl(Ax.XTick);
Ax.XTickLabelRotation = 45;
xlabel("Time")
ylabel("Value")
title("T1: Boxchart")
axis("padded")
% get(hb)
Adi Purwandana
on 19 Oct 2024
Edited: Adi Purwandana
on 19 Oct 2024
Also, if possible... (I tried but still got no solution), can we put a mean value plot connecting the boxplot? There's likely I need groupsummary.
Just like:
Thank you and sorry for so many following questions.
Everything in my previous code should have worked in R2022a. I checked — all those functions existed from at least R2020b, most earlier. Apparently boxchart wasn’t compatible with string arrays then, so we’ll see if a cell array works with it this time. Adding the mean line to the boxchart is straightforward. I added a legend as well to both plots. You can reposition each as desired, or not use it if you don’t want it.
Try this —
T1 = readtable('datamine.xlsx')
idx = (T1.value_y >= 12) & (T1.value_y < 13.4); % selecting only y-value with a certain range
T1 = T1(idx,:)
TT1 = table2timetable(T1);
SEM = @(x) std(x)/sqrt(numel(x));
TT1momean = retime(TT1, 'monthly', 'mean');
TT1mosem = retime(TT1, 'monthly', SEM);
TT1monum = retime(TT1, 'monthly', 'count');
[Nmin,Nmax] = bounds(TT1monum{:,1})
% TT1Time = TT1mosem.time;
% TT1SEM = TT1mosem{:,1};
figure
hplt = plot(TT1momean.time, TT1momean{:,1}, '-k', "DisplayName",'Monthly Mean');
hold on
hpch = patch([TT1mosem.time; flip(TT1mosem.time)], [TT1momean{:,1}-TT1mosem{:,1}*1.96; flip(TT1momean{:,1}+TT1mosem{:,1}*1.96)], 'r', 'FaceAlpha',0.375, 'EdgeColor','r', 'EdgeAlpha',0.5, "DisplayName","±95% Confidence Interval");
hold off
grid
xlabel('Time')
ylabel('Value')
title('Value Mean & ±95% CI')
legend("Location","NE")
[GV, ID1, ID2] = findgroups(year(T1.time), month(T1.time)); % Grouping Variable For ‘boxchart’
% [b1,b2] = bounds(ID1)
xtl = compose('%s',datetime(ID1,ID2,ones(size(ID1)), Format="MMM yyyy")); % Cell Array Of Time Values
% disp(xtl)
figure
hb = boxchart(GV, T1.value_y, "DisplayName","Boxchart of Monthly Data");
hold on
plot((1:height(TT1momean)), TT1momean{:,1}, '.-r', "DisplayName","Monthly Mean")
hold off
hb.JitterOutliers = "on";
hb.MarkerStyle = ".";
hb.Notch = "on";
hb.BoxMedianLineColor = "r";
Ax = gca;
xt = Ax.XTick;
Ax.XTick = [1 xt];
Ax.XTickLabel = xtl(Ax.XTick);
Ax.XTickLabelRotation = 45;
xlabel("Time")
ylabel("Value")
title("Boxchart of Monthly Data")
hlgnd = legend("Location","northoutside", "Orientation","horizontal");
axis("padded")
% get(hb)
That should work.
.
Adi Purwandana
on 20 Oct 2024
I really appreciate your help. FYI, I just installed 2024b and it works fine, but not for the 2022a.
It is mostlikely those (adjustable BoxMedianLineColor and date format syntax) are not available in 2022a. I don't know. Anyway, thank you for such very helpful solution. The best as always!
Star Strider
on 20 Oct 2024
As always, my pleasure!
I did not consider that the other options for boxchart could have been the problem. The ‘Version History’ in the boxchart documentation did not mention any changes.
Adi Purwandana
on 21 Oct 2024
Yes.
Create a vector using the colon, : operator, with the second (‘increment’ or ‘step’) value being whatever integer value you want. Here, it needs to be an integer since it will be used as a subscript, however that is the only constraint:
Ax.XTick = [1 : 3 : xt(end)]; % Every 3 Months
Ax.XTick = [1 : 6 : xt(end)]; % Every 6 Months
In this code, I commented the every 6 month line so it would not run.
It can also be whatever months you want to display, for example:
Ax.XTick = [2 3 4 7 10 17 22 28 33]; % Specific Months
The only constraint is that they be integers from 1 to 35 (in this instance).
Try this —
T1 = readtable('datamine.xlsx')
idx = (T1.value_y >= 12) & (T1.value_y < 13.4); % selecting only y-value with a certain range
T1 = T1(idx,:)
TT1 = table2timetable(T1);
SEM = @(x) std(x)/sqrt(numel(x));
TT1momean = retime(TT1, 'monthly', 'mean');
TT1mosem = retime(TT1, 'monthly', SEM);
TT1monum = retime(TT1, 'monthly', 'count');
[Nmin,Nmax] = bounds(TT1monum{:,1})
% TT1Time = TT1mosem.time;
% TT1SEM = TT1mosem{:,1};
figure
hplt = plot(TT1momean.time, TT1momean{:,1}, '-k', "DisplayName",'Monthly Mean');
hold on
hpch = patch([TT1mosem.time; flip(TT1mosem.time)], [TT1momean{:,1}-TT1mosem{:,1}*1.96; flip(TT1momean{:,1}+TT1mosem{:,1}*1.96)], 'r', 'FaceAlpha',0.375, 'EdgeColor','r', 'EdgeAlpha',0.5, "DisplayName","±95% Confidence Interval");
hold off
grid
xlabel('Time')
ylabel('Value')
title('Value Mean & ±95% CI')
legend("Location","NE")
[GV, ID1, ID2] = findgroups(year(T1.time), month(T1.time)); % Grouping Variable For ‘boxchart’
% [b1,b2] = bounds(ID1)
xtl = compose('%s',datetime(ID1,ID2,ones(size(ID1)), Format="MMM yyyy")); % Cell Array Of Time Values
% disp(xtl)
figure
hb = boxchart(GV, T1.value_y, "DisplayName","Boxchart of Monthly Data");
hold on
plot((1:height(TT1momean)), TT1momean{:,1}, '.-r', "DisplayName","Monthly Mean")
hold off
hb.JitterOutliers = "on";
hb.MarkerStyle = ".";
hb.Notch = "on";
hb.BoxMedianLineColor = "r";
Ax = gca;
xt = Ax.XTick;
Ax.XTick = [1 : 3 : xt(end)]; % Every 3 Months
% Ax.XTick = [1 : 6 : xt(end)]; % Every 6 Months
Ax.XTickLabel = xtl(Ax.XTick);
Ax.XTickLabelRotation = 45;
xlabel("Time")
ylabel("Value")
title("Boxchart of Monthly Data")
hlgnd = legend("Location","northoutside", "Orientation","horizontal");
axis("padded")
% get(hb)
.
Adi Purwandana
on 21 Oct 2024
Star Strider
on 21 Oct 2024
As always, my pleasure!
More Answers (0)
Categories
Find more on Time Series Events in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
