Multilevel logitic regression power test

[Bianca Onita]

Hello,

I'm having difficulty finding materials on calculating test power (analysis power) for multilevel logistic regression analysis. I specifically need to calculate the power for the fixed coefficients. I noticed that STATA itself has commands for calculating sample size and test power. Does anyone know if there is a command for calculating test power for multilevel logistic regression?

[Joseph Coveney]

You might need to resort to simulation.

You don't mention any particulars, and so I cannot suggest anything specific.

But if you Google power simulation site:stata.com, it'll bring up a number of hits on the general technique, including an FAQ on the company's website, a series of how-to articles on its blog site and references to the custom (user) method of its official Stata power command.

[Bianca Onita]

Thanks for the reply Joshep! I'm going to try your suggestion.

In my case, I need to calculate the power a posteriori because the sample size was already determined. Some informations about my situation are: I have a logistic regression multilevel model, with 2 levels. Level 1 has n=2842 and level 2 has n=1421. I have multiple independent variables in my model.

[Joseph Coveney]

I need to calculate the power a posteriori because the sample size was already determined.

You mentioned the counts at each level, but power or sample-size estimation requires specifying a null- and alternative-hypothesis pair, too. In addition, if you're going to use simulation, you'll need to consider the variance at the higher level. You repeat the simulation exercise across a range of plausible values of the random effects variance.

[Bianca Onita]

Thank you, Joseph!

I was trying to run the simulations, but the command is not running in my STATA and I don't know why. As especified in https://blog.stata.com I wrote a program to run my simulations with power in the model melogit, but I always get the message:
power: invalid method simmelogit
The simmelogit method is not one of the officially supported power methods.

Does anyone know why?

[Joseph Coveney]

. . . I wrote a program to run my simulations with power . . .

I don't know why your program doesn't work with power, but that's not how I would approach this.

Instead I would write a simulation program that takes as arguments those parameters that you wish to explore, for example, effect size, sample size (number of higher-level entities) and cluster size (number of lower-level observations per entity) and as I mentioned above magnitude of the variance at the higher level. I show an example of this in the code section below. For illustration I've chosen a minimum detected effect size of an odds ratio of 2 (in your case with multiple predictors it would be an adjusted odds ratio of 2) centered about odds of 1:1 (50% outcome rate), and sample size throughout a range with a fixed cluster size of 3 (which can be varied as an option) and a range of values of the random effects variance expressed in terms of the corresponding intraclass correlation coefficient (ICC).

version 18.0

clear *

// seedem
set seed 1172256435

program define or2, rclass
    version 18.0
    args x

    return scalar fx = (`x' / (1 - `x') / ((1 - `x') / `x') - 2)^2
end

minbound or2, range(0.55 0.60)
tempname p0 p1
scalar define `p0' = 1 - r(x)
scalar define `p1' = r(x)

program define simEm, rclass
    version 18.0
    syntax , p0(name) p1(name) [n(integer 250) k(integer 3) icc(name)]

    drop _all
    quietly set obs `n'
    generate int pid = _n
    generate double pid_u = rnormal(0, sqrt(`icc' * _pi^2 / (3 - 3 * `icc')))

    generate byte trt = mod(_n, 2)

    expand `k'
    generate double xbu = logit(cond(trt, `p1', `p0')) + pid_u
    generate byte out = rbinomial(1, invlogit(xbu))

    xtlogit out i.trt, i(pid) re
    return scalar pos = r(table)["pvalue", "out:1.trt"] < 0.05
end

frame create PowerAnalysis double icc int n double power

tempname icc
foreach icc100 in 20 50 80 {

    scalar define `icc' = `icc100' / 100
    display in text _newline(1) "ICC = " %04.2f scalar(`icc')

    forvalues n = 200(100)400 {
        quietly simulate pos = r(pos), reps(100): simEm , p0(`p0') p1(`p1') ///
            n(`n') icc(`icc')
        summarize pos, meanonly
        display in smcl as text "N = " as result "`n'" ///
            as text " power = " as result %04.2f r(mean)
        frame post PowerAnalysis (`icc') (`n') (r(mean))
    }
}

cwf PowerAnalysis
local opts sort lcolor(black) mcolor(black) mfcolor(white)
graph twoway ///
    connected power n if icc == 0.2, `opts' msymbol(O) msize(medium) || ///
    connected power n if icc == 0.5, `opts' msymbol(T) msize(medlarge) || ///
    connected power n if icc == 0.8, `opts' msymbol(S) msize(medium) ///
        xtitle(Sample Size) xlabel(200(100)400) ///
        ytitle(Power) ylabel(0(0.2)1, format(%04.2f) angle(horizontal) nogrid) ///
        scheme(s2color) legend(off)
graph export PowerAnalysis.png

exit

You can see (output listed below and displayed in the attached graph) that power strongly depends upon magnitude of the variance component, and I'm not sure that that can be easily specified via the official power command.:

ICC = 0.20
N = 200 power = 0.89
N = 300 power = 0.97
N = 400 power = 0.98

ICC = 0.50
N = 200 power = 0.57
N = 300 power = 0.73
N = 400 power = 0.80

ICC = 0.80
N = 200 power = 0.23
N = 300 power = 0.39
N = 400 power = 0.48

I've attached the do-file and log-file if you're interested in pursuing this approach further with your study's specifics.

[Bianca Onita]

Thank you so much for your detailed explanation and your attention to my difficulties!

I am going to try these lines of commands that you kindly attached and see if I succeed in my attempts!

Thank you again, Joseph!