Obtaining overall slope and level 2 slope residual from multilevel model (xtmixed)

[Anika Knueppel]

Hi, I need help with ensuring I'm using the right postestimations commands to obtain a slope through repeated measurements of a biomarker using a multilevel model.
The data structure is in a longform with per study participant 3 biomarker measurements over the day with exact time of measurement for each.
I am following a model used in a research paper described which can be paraphrased as: The slope across the day was derived from a hierarchical linear model (multilevel model) to predict the biomarker across the measurements with measurement occasion used as level 1 identifier, participant as level 2 identifier, sample time as independent variable and random intercepts and random slopes. The slope was estimated per participant as the overall slope + level 2 slope residual.

id=participant identifier
occasion=measurement occasion 1, 2 or 3 (based on the three measurements taken over the day
time= exact time of day of measurement
biomarker= value of biomarker measurement

I have been running the following:
xtset id occasion
xtmixed biomarker time || id:, mle cov(unstr) variance
predict level2s, reffects level(id)
gen slope=_b[measurement-time]+level2s

It would be great to hear if this looks correct to anyone?
Thank you so much in advance

[Joseph Coveney]

A couple of comments.

First, I don't know what version of Stata you're using, but the name of the command currently is mixed. It has maximum likelihood as the default method and displays the variance components' variance by default and so neither of those options need to be specified. (I assume that your sample size justifies using that method over, say, a combination of REML and Kenward-Roger options.)

Second, xtset doesn't do anything in this context, and you can omit that line.

Third, you haven't included your time variable in your random effects equation, and so there won't be any random slope fitted. When you do include the time variable in the random effects equation, you'll need to specify the prediction variables individually or as a stub*. (By the way, just what is your time variable's name? You have it as "occasion" in the first line of code and then "time" in the second.)

So, maybe your snippet of code would be better as something more like the following.

/* xtset id occasion */
mixed biomarker c.time || id: time, cov(unstructured)
predict double blups*, reffects // level(id)
generate double slope _hat = _b[time] + blup1

I assume that you know about shrinkage. Run the following code and see how close the random slope values are to the actual in a set up (three observations per participant) similar to yours. (You can adjust the magnitudes of the fixed effect slope and covariance matrix of the random parameters by the observed values from your fitted model in order to more closely mimic your situation.)

version 17.0

clear *

// seedem
set seed 1484975132

// Participants
quietly set obs 250
generate int pid = _n

// Random intercept and slope
drawnorm i s, double corr(1 -0.25 \ -0.25 1)

// Occasions
quietly expand 3
bysort pid: generate byte tim = _n

// Biomarker
generate double bmk = ///
    10 + i + ///
     (tim - 1.5) / 1.5 * (-1 + s) + ///
        rnormal(0, sqrt(0.25))

*
* Begin here
*
mixed bmk c.tim || pid: tim, covariance(unstructured) nolrtest nolog
predict double blups*, reffects

foreach var of varlist blups? {
    local varlabel : variable label `var'
    display in smcl as text "`var': `varlabel'"
}

// Shrinkage
graph twoway ///
    line s s, sort lcolor(black) lpattern(dash) || ///
    scatter blups1 s if tim == 1, mcolor(black) msize(vsmall) ///
        ylabel( , angle(horizontal) nogrid) legend(off)

// Try this
generate double s_hat = _b[c.tim] + blups1

exit