Help with using Monte Carlo Simulation to do sensitivity power analysis on a mixed multilevel model

[Apostolos Ng]

Hi,
I am a psychology honours student who is using a mixed multilevel model for one of my research questions. I am considering using Monte Carlo simulation to do a sensitivity power analysis. However, we have not learn this at uni, and I am not sure if I am doing it correctly. Would someone be able to check if my coding is correct? It seems to run fine on my version of Stata 17, but as I am unfamiliar with it I just wanted to check if I have done the coding correctly. My variables are
pid: participant ID (level 2)
agegrp: age group, 1 = older, 2 = younger (older is reference level)
cval: cue valence, 1 = negative, 2 = neutral, 3 = positive
avg_rating_cue: affective appraisal rating (DV), averaged within each cue valence condition
prop_spec_val_c: proportion of specific memories in each valence condition (predictor, grand centred)

Step 1: Estimate inputs from the real data
* SD components of the predictor (valence main effect removed)
mixed prop_spec_val_c i.cval || pid:, reml
scalar sd_spec_b = exp(_b[lns1_1_1:_cons]) // between-person SD
scalar sd_spec_w = exp(_b[lnsig_e:_cons]) // within-person SD
scalar sd_spec_tot = sqrt(sd_spec_b^2 + sd_spec_w^2)

* Outcome model (observed)
mixed avg_rating_cue i.agegrp##c.prop_spec_val_c i.cval || pid:, reml
scalar sd_intercept = exp(_b[lns1_1_1:_cons]) // random-intercept SD
scalar sd_resid = exp(_b[lnsig_e:_cons]) // residual SD
scalar b_cons = _b[_cons]
scalar b_agegrp = _b[2.agegrp]
scalar b_spec = _b[c.prop_spec_val_c] // slope in the reference group (older)

Step 2: simulate one dataset and test the interaction
capture program drop sim_rq3
program sim_rq3, rclass
version 16
syntax, delta(real) n1(integer) n2(integer) nobs(integer) [alpha(real 0.05)]
local ntot = `n1' + `n2'
local reject = .
local conv = 0
quietly {
drop _all
set obs `ntot'
generate long pid = _n
generate byte agegrp = cond(_n <= `n1', 1, 2) // 1 = Older, 2 = Younger

generate double u_i = rnormal(0, sd_intercept) // random intercept
generate double s_i = rnormal(0, sd_spec_b) // person's mean specificity
expand `nobs' // 3 observations per person
generate double spec = s_i + rnormal(0, sd_spec_w)

* delta = true slope difference (added to group 1)
generate double slope = cond(agegrp==1, b_spec + `delta', b_spec)
generate double rating = b_cons + b_agegrp*(agegrp==2) ///
+ slope*spec + u_i + rnormal(0, sd_resid)

capture mixed rating i.agegrp##c.spec || pid:, reml iter(200)
if _rc == 0 {
if e(converged) == 1 {
local conv = 1
testparm i.agegrp#c.spec
local p = chi2tail(r(df), r(chi2)) // Wald test
local reject = (`p' < `alpha')
}
}
}
return scalar reject = `reject'
return scalar conv = `conv'
end

Step 3: power curve
local n_older = 34
local n_younger = 36
local nobs = 3
local reps = 1000
local seed = 12345

tempname results
postfile `results' double delta double delta_per_sd double power double conv_rate ///
using rq3_power_curve.dta, replace

foreach d of numlist 0 1.0(0.1)1.4 {
quietly simulate reject=r(reject) conv=r(conv), reps(`reps') seed(`seed') nodots: ///
sim_rq3, delta(`d') n1(`n_older') n2(`n_younger') nobs(`nobs')
quietly summarize reject
local pwr = r(mean)
quietly summarize conv
local cr = r(mean)
post `results' (`d') (`d'*sd_spec_tot) (`pwr') (`cr')
}
postclose `results'

Thank you for your help.

[Joseph Coveney]

I am considering using Monte Carlo simulation to do a sensitivity power analysis. However . . . I am not sure if I am doing it correctly. Would someone be able to check if my coding is correct?

While waiting for someone to jump in who's more able to follow what you're trying to do, I'll give it a shot with a couple of questions for clarification.

1. The purpose of this is "a sensitivity power analysis". How does that differ from a widely disparaged post-hoc power analysis?

2. The hierarchical regression model's outcome (response) variable, "affective appraisal rating", is "averaged within each cue valence condition". So, each study participant is observed multiple times under each of the three "cue valence conditions" ("negative", "neutral" and "positive")? If so, then wouldn't you want to have that at an additional level in your hierarchical regression model in order to capture that contribution to the uncertainty estimates?

3. You haven't really described your study's design and haven't shown any example of the dataset, but if I've got the gist of your study correct, then your hierarchical regression model's interaction term (i.agegrp##c.prop_spec_val_c) includes both a participant-invariant predictor (age category) and a predictor that varies within each participant ("proportion of specific memories in each valence condition" over three "cue valence conditions" for each participant), which represents a so-called cross-level interaction. If that's the case, then the widespread recommendation is to include a random slope for the participant-varying predictor that's a component of the cross-level interaction term.

4. The (fixed) sample size I guess is determined by what your study happens to have enrolled in each age category, and their sum is only 70 study participants. If that's the case, you might want to add the option to your hierarchical regression model for small-sample degrees of freedom adjustment, say, dfmethod(kroger).

[Joseph Coveney]

I forgot to include one:

5. I don't understand the point of your Step 1, but it implies that "proportion of specific memories" is not exogenous, that is, the variable that you're using as a predictor in your hierarchical regression model is itself an outcome (response) variable, perhaps an intermediate one. Is that correct?