Data visualisation and mixed model (individual) effect estimates for plotting?

[patrick handcock]

Hi Statalisters,

I was hoping someone might be able to help with this. I have run a mixed model analysis of a cross-over experimental study looking at the effect of activity breaks vs. sitting on some acute biomarker responses (GLP1 and GIP iAUC's). I have presented the overall response to each of the three conditions, adjusted for a few covariates and baseline variables, in bar charts (or 'detonator plots'). A reviewer has suggested I used a different plot to better represent the data (I guess by showing individual data points so as to be more transparent (e.g. they referred to this and this paper).

I agree this would be a good thing to do and a more transparent way to present the data, but I am wondering how best to do this on the modelled data (i.e. with models adjusting for other variables like age, bmi, etc). I was thinking I could show the raw individual iAUC values for each participant per condition as dot points on the plot to show spread of data, then still also present the mean effect and SEM (or 95% CI), but that doesn't seem quite right given the mean effect is based on modelled data (adjusting for covariates and baseline values, etc). I am wondering if this modelled data can be pulled out from the mixed model, or best just to find a different way to present the individual data and not just the detonator plot?

CONTEXT
Stats section:
Postprandial GLP-1 and GIP responses during each trial were summarised as mean values and total and net incremental area under the curve (tAUC, iAUCnet), calculated using the trapezoidal method. Generalised linear mixed-models with random intercepts examined the differential effects of the experimental conditions on all summary outcomes using Stata 14 (Stata-Corp LP), adjusting for potential covariates explaining residual outcome variance (age, BMI, sex), baseline values, and period effects (treatment order). Sex-by-condition interactions were tested for each outcome.

Results
Figure 1 shows mean plasma GLP-1 and GIP concentrations across the three conditions (n=23 participants; 14 males; 9 females). Overall (6.5 hours) and meal-specific tAUC and iAUCnet were significantly lower in LW and SRA versus SIT for GIP, but not for GLP-1 (Figure 1; Table 1). No significant differences were observed between the LW and SRA conditions.

Mixed model Example CODE to compare the conditions. Is there a way to get individual data estimates to also plot perhaps behind the bar charts?....

*******************venous GIP analysis**************************
mixed ALLtotalAUC_gip_hr i.condition gip_baseline age i.sex bmi i.trialorder || id: , mle var
margins i.condition
pwcompare i.condition, mcompare(noadjust) effects


** GENDER INTERACTION **
mixed ALLtotalAUC_gip_hr i.condition##i.sex gip_baseline age bmi i.trialorder || id: , mle var
testparm i.condition##i.sex
margins i.condition##i.sex
pwcompare i.condition##i.sex, mcompare(noadjust) effects


Thanks!

[Erik Ruzek]

A couple of thoughts, both about your models and then about your graphing question.

Regarding the models, mixed models for continuous outcomes rely on large sample theory for both the variance of the random effects and the standard errors of the fixed effects. With only 23 participants, you should be using restricted maximum likelihood estimation (mixed option reml), ideally with Kenward-Roger standard error corrections (mixed option, dfmethod(kroger)). See here and also McNeish's (2017) article on the topic of REML.

One more model-based question. There appears to be a clear upward trend over time in plot C. That suggests that you may want to include a variable that indexes the passage of time (however it is most relevant in your study) as a predictor. To be transparent, I know nothing about your outcome, I just find it strange that you wouldn't adjust for time effects given such a plot.

With regards to your graphing issue, I wonder if the reviewer understands that the detonator plot is not an unadjusted plot of the means? Is it clear in the in-text description of the plot that the estimates are based on predictions from your model? That is not obvious in the Figure text. You can certainly overlay data points onto a marginsplot.

use https://www.stata-press.com/data/mlmus3/faculty
sum market
gen C_market = market - r(mean)
regress salary C_market i.male
margins male, at(C_market = (-.2(.2).4))
marginsplot, addplot(scatter salary C_market if male==0, msymb(Oh) mcolor(gs2) || scatter salary C_market if male==1, msymb(O) mcolor(gs10) color(%40))

[patrick handcock]

Thanks so much for your astute and helpful thoughts Erik Ruzek .

1. RE. using REML. Good point - I will make this update to the model - although if I recall correctly it did not impact the results when I checked (can do anyway to be conservative).

2. Yes using a time variable on the model is another way to look at this (and perhaps more robust), but the iAUC variable already captures/summarises the 'incremental or total area under the curve' (ie. the meal-specific response of the analyte) over time, so in effect does a similar thing and works better with such outcomes. It is quite commonly used in acute experimental studies - perhaps as it is more sensitive to picking up dynamic differences in 'response' over time.

3. Thanks for this RE. detonator plot. Yes my initial thought was to respond and make it clearer that the detonator plot is indeed adjusted/modelled (inference), but I wondered if even so whether that precluded me from also adding the individual data points for transparency on spread, etc. What I was less clear on was whether I should overlay the raw unadjusted iAUC data, or try and somehow present the adjusted data from that model?... I am doing the graphing in GraphPad Prism - do you know if there is a way to pull out the adjusted individual data points to overlay?

4. On point 3 - I also wondered if perhaps presenting a box plot of some other kind of plot that shows spread might be even better (e.g. violin or indeed 'spread' plots) - but then I wondered if that fundamentally misses the point that this is modelled data - so showing the spread is less relevant (as modelled data is about inference)?... Still, I think the transparency still helps so spread might be most ideal - but mixing raw individual data and modelled mean (95% CI) did not seem quite right?... (sorry repeating point 3 a bit here)

Thanks again for your time!

[Tiago Pereira]

Erik's suggestions are incredibly important.
--
Regarding the plot: You can use spaghetti plots, showing the observed values. Alternatively, you can plot the predicted values from -mixed-, encompassing the linear prediction plus + predicted random effect.

mixed ALLtotalAUC_gip_hr i.condition##i.sex gip_baseline age bmi i.trialorder || id: , reml    dfmethod(kroger) 

predict y_hat, fitted

See these posts for further insights:
https://www.statalist.org/forums/for...oup-mean-lines
https://errickson.net/stata.html

See the manual for mixed for further information on the predictions:

https://www.stata.com/manuals13/memi...estimation.pdf

Nice example from the internet showing individual trajectories and the overall trend (probably not created in Stata, but it would be easy to implement it in Stata). I am not telling you this is the best way to present your data, but it could be one of the possible strategies.

Source: https://www.storytellingwithdata.com...paghetti-graph

[patrick handcock]

Thanks a lot for all the great advice! Much appreciated!