Twoway graph with 2 y- axes on time series data

[Neshika Samarasekera]

My data comes from an observational study which measured cell counts (lymc) and brain swelling (total_pho) in 71 patients at the same time points after stroke. It measured both variables at the same 4 time points (numbered 1,2,3,4) after the stroke. I am trying to plot the 2 variables on the same graph -both variables should be line graphs and I would like the median of each variable plotted at a given time, but the y axes will have different scales. The data is in a long format (sample posted below) - and there are some missing values e.g. person 134 doesn't have cell count or brain swelling measured at times 3 and 4.

id time total_pho lymc
131 1 2.13 1.32
131 2 3.68 2.08
131 3 2.41 2.21
131 4 2.21
132 1 59.96 1.31
132 2 70.5 2.58
132 3 98.08 2.07
132 4 114.17 1.99
133 1 24.06 3.1
134 2 29.1 2.27
134 3
134 4
135 1 51.8 .79
135 2 50.43 1.11

I used this code -

graph twoway (line lymc time, yaxis(1)) (line total_pho time, yaxis(2))

It produced the graph with the correct axes but doesn't seem to recognise that the data is time series, so has produced a separate line for each separate participant.

How do I sort that? Is there also a way to plot the interquartile range for each varible at each of the time points?

thankyou

[Andrew Musau]

You have panel data. If you want to plot the average trends, you'd need to collapse the data first. Color coding the axes titles would help with interpretation in your design.

clear
input float(id time total_pho lymc)
131 1 2.13 1.32
131 2 3.68 2.08
131 3 2.41 2.21
131 4 2.21 .
132 1 59.96 1.31
132 2 70.5 2.58
132 3 98.08 2.07
132 4 114.17 1.99
133 1 24.06 3.1
134 2 29.1 2.27
134 3 . .
134 4 . .
135 1 51.8 .79
135 2 50.43 1.11
end

preserve
collapse lymc total_pho, by(time)
graph twoway (line lymc time, yaxis(1) ytitle(, color(stc1))) ///
(line total_pho time, yaxis(2) ytitle(, color(stc2) axis(2)) leg(off))
restore

[Andrew Musau]

See

help collapse

for more statistics other than the mean.

[Nick Cox]

This is similar in spirit to Andrew Musau's suggestion, but different in detail.

I note that you want to show IQR, which I take it is best done implicitly by showing quartiles as well as median. So, that implies the box of a box plot, but nothing beyond. With 71 patients you can but possibly also should show much more detail, but leave that on one side

Missing values also raise the possibility that only patients with complete profiles should be plotted.

People often ask for two (or even more) different y axes in the same plot. I am sceptical that such displays are easy to understand, let alone effective, especially you are showing three statistics for each time and measure.

Here's a different take.

clear 
input id time total_pho lymc
131 1 2.13 1.32
131 2 3.68 2.08
131 3 2.41 2.21
131 4 2.21
132 1 59.96 1.31
132 2 70.5 2.58
132 3 98.08 2.07
132 4 114.17 1.99
133 1 24.06 3.1
134 2 29.1 2.27
134 3
134 4
135 1 51.8 .79
135 2 50.43 1.11
end 

preserve 
drop id 

capture set scheme stcolor 

stack time total_pho time lymc, into(time value) clear 
label def _stack 1 "brain swelling" 2 "cell counts"
label val _stack _stack 

foreach p in 25 50 75 { 
egen p`p' = pctile(value), p(`p') by(time _stack)
} 

egen n = count(value), by(time _stack)
gen toshow = "{it:n} = " + strofreal(n)

egen min = min(value), by(_stack)
egen max = max(value), by(_stack)
gen where = min - (max - min) * 0.1 

scatter p50 time , ms(Dh) by(_stack, col(1) note("") yrescale legend(off)) ///
subtitle(, pos(9) nobexpand nobox fcolor(none )) ///
|| rbar p25 p75 time, barw(0.2) fcolor(none) /// 
|| scatter where time, ms(none) mlabel(toshow) mlabsize(large) mlabpos(0) mlabc(black)

restore

[Neshika Samarasekera]

thank you so much for taking the time to respond. I will explore both of these options

[Neshika Samarasekera]

Dear Andrew, Nick

Andrew - I found your suggestion particularly helpful, and used the following code which worked (Graph 1, collapsed data also attached) - thank you ++ Nick I also tried yours out so thank you for taking the time.

graph twoway (line medlymc time, yaxis(1) ytitle(median lymphocyte count, color(red))) (line medtotphotime, yaxis(2) ytitle(median PHO volume, color(blue) axis(2)) leg(off))

I have been trying different ways of changing the graph in two ways - 1 - I would like to have the interquartile range at each time point for both brain swelling (PHO volume) and cell count (lymphocyte count). I have been using the command rarea but can't get this to work, as stata seems to confuse the graph so that the axis for lymphocyte count changes from 0-60 and the tendency of both variables to increase isn't clear (Graph 2) I also wanted to specify the scale for the y axes so the lymphocyte count goes from 0-2 with intervals of 0.5 and PHO volume (brain swelling) goes from 0-70 with intervals of 10. How do I do that? I tried ylabel but couldn't do it.
These are the codes I used

graph twoway (line medlymc time, yaxis(1) ytitle(median lymphocyte count, color(red))) (line medtotpho time, yaxis(2) ytitle(median PHO volume, color(blue) axis(2)) leg(off)) rarea p75lymc p25lymc time||rarea p75totpho p25totpho time||

graph twoway (line medlymc time, yaxis(1) ylabel(0(0.5)2) ytitle(median lymphocytecount, color(red))) (line medtotpho time, yaxis(2) ylabel(0(10)70) ytitle(medianPHO volume, color(blue) axis(2)) leg(off)) rarea p75lymc p25lymc time||rarea p75totpho p25totpho time||


Thank you for your time, I do appreciate this

Nesh

[Nick Cox]

Please don't post .gph or .dta attachments, as explained in the FAQ Advice. https://www.statalist.org/forums/help#stata

Please post images as ,png attachments and data examples as code using dataex.

[Nick Cox]

Here's your data. From a sight of your graphs, you're not using Stata 19. You should tell us what version you're using if not.

* Example generated by -dataex-. For more info, type help dataex
clear
input byte time float(medtotpho medlymc p25totpho p25lymc p75totpho p75lymc)
1  10.01  1.32  6.85  .92  32.59 1.86
2 18.425 1.325  9.57 1.05     43 1.75
3  20.11  1.44 10.01 1.05 49.777 1.81
4  26.01  1.64  11.9 1.14  62.07    2
end
label var medtotpho "(p 50) total_pho" 
label var medlymc "(p 50) lymc" 
label var p25totpho "(p 25) total_pho" 
label var p25lymc "(p 25) lymc" 
label var p75totpho "(p 75) total_pho" 
label var p75lymc "(p 75) lymc"

I don't think Stata is confused by your graph code. It's just giving you what you asked for. You need with your design to tell graph that median and quartiles for one variable are all on the same scale and those for the other variable are all on a different scale.

Here is first the graph I still think works better for you. Then I worked like this. The pho variables belong in a spece from 5 to 65 (say) and the lymc variables in a space from 0.8 to 1.2. So we scale to the interval [0, 1]. Then we fix the axis labels using mylabels from the Stata Journal.

* Example generated by -dataex-. For more info, type help dataex
clear
input byte time float(medtotpho medlymc p25totpho p25lymc p75totpho p75lymc)
1  10.01  1.32  6.85  .92  32.59 1.86
2 18.425 1.325  9.57 1.05     43 1.75
3  20.11  1.44 10.01 1.05 49.777 1.81
4  26.01  1.64  11.9 1.14  62.07    2
end
label var medtotpho "(p 50) total_pho" 
label var medlymc "(p 50) lymc" 
label var p25totpho "(p 25) total_pho" 
label var p25lymc "(p 25) lymc" 
label var p75totpho "(p 75) total_pho" 
label var p75lymc "(p 75) lymc" 

twoway rarea p75lymc p25lymc time, lcolor(red) fcolor(red%10) ///
|| connected medlymc time, ytitle(lymphocyte count, color(red)) ///
xsc(off) fysize(45) legend(off) name(G1, replace)

twoway rarea p75totpho p25totpho time, fysize(55) lcolor(blue) fcolor(blue%19) ///
|| connected medtotpho time, ytitle(PHO volume, color(blue)) color(blue)  ///
legend(off) name(G2, replace)

graph combine G1 G2, xcommon col(1) name(G3, replace) 

foreach v of var *lymc { 
    gen `v'2  = (`v' - 0.8)/1.2
}

foreach v of var *pho {
    gen `v'2 = (`v' - 5)/60
}

su *2

mylabels 10(10)60, myscale((@ - 5)/60) local(right)
mylabels 0.8(0.2)2, myscale((@ - 0.8)/1.2) local(left)

twoway rarea p75lymc2 p25lymc2 time, lcolor(red) fcolor(red%10) ///
|| rarea p75totpho2 p25totpho2 time, lcolor(blue) fcolor(blue%10) ///
|| connected medlymc2 time, yaxis(1 2) color(red) ytitle(lymphocyte count, color(red) axis(1)) yla(`left', labcolor(red) axis(1)) ///
|| connected medtotpho2 time, color(blue) ytitle(PHO volume, color(blue) axis(2)) color(blue)  ///
yla(`right', labcolor(blue) axis(2)) legend(off) name(G4, replace)

[Nick Cox]

Extra code in the same spirit for your specified axis labels:

foreach v of var *lymc { 
    replace `v'2  = `v'/2
}

foreach v of var *pho {
    replace `v'2 = `v'/70
}

mylabels 0(10)70, myscale(@/70) local(right)
mylabels 0(0.5)2, myscale(@/2) local(left)

twoway rarea p75lymc2 p25lymc2 time, lcolor(red) fcolor(red%10) ///
|| rarea p75totpho2 p25totpho2 time, lcolor(blue) fcolor(blue%10) ///
|| connected medlymc2 time, yaxis(1 2) color(red) ytitle(lymphocyte count, color(red) axis(1)) yla(`left', labcolor(red) axis(1)) ///
|| connected medtotpho2 time, color(blue) ytitle(PHO volume, color(blue) axis(2)) color(blue)  ///
yla(`right', labcolor(blue) axis(2)) legend(off) name(G5, replace)

[Nick Cox]

Naturally, you should want more explanation of units of measurement. In my view, a legend would be over the top, as well as taking up far too much space. I would recommend putting explanations of median and quartiles as extra titles (or note or caption).

[Neshika Samarasekera]

thank you Nick++

I am using Stata 11.2 (rather an old version I know) Does that make a difference to the code I am running? presume it does. I have run all the updates on my version of it

Your first graph looks like exactly what I need- I will try to model that one. thank you, and apologies for posting those attachments

[Nick Cox]

Telling us the version you're using (if it is not the latest) is a longstanding request in the FAQ Advice all posters are asked to read before posting.

https://www.statalist.org/forums/help#version

I don't have access to Stata 11 any more, but I am optimistic that the code in #8 and #9 should work in 11. You would need to install mylabels but if that is impossible work-arounds can be imagined.

In Stata 11, I strongly recommend using almost any scheme but the default s2color.

set scheme s1color

is one possibility.

 search mylabels, sj

Search of official help files, FAQs, Examples, and Stata Journals

SJ-24-1 gr0092_1  . . . . . . . . . . . . . . . . Software update for mylabels
        (help nicelabels, mylabels, myticks if installed) . . . . .  N. J. Cox
        Q1/24   SJ 24(1):182--184
        fixes a bug that could bite if the options myscale() and
        clean were specified together

SJ-22-4 gr0092  . . . . . . . . . . . . Speaking Stata: Automating axis labels
        (help nicelabels, mylabels, myticks if installed) . . . . .  N. J. Cox
        Q4/22   SJ 22(4):975--995
        provides commands to handle two common problems with graph
        axis labels: decide in advance on some "nice" numbers to
        use on one or both axes and show particular labels on some
        transformed scale

[Neshika Samarasekera]

thank you Nick.