Cook's D after meglm

[Amar Singh]

Hello all,

I wanted to generate Cook's D after meglm command, so I used I used 'predict cooksd, cooksd', but it says the option is not allowed. Do you think there is anyway we can calculate Cook's D after meglm? Please help. Thank you.

[Joseph Coveney]

Well, here's a brute-force way for the fixed effects equation. The example uses a logistic regression model with a single random effect for illustration, but you can modify it for other distribution families, link functions and random effects equations.

version 14.1

clear *
set more off
set seed 1349773

quietly set obs 250
generate int s = _n
generate double u = rnormal()

forvalues t = 1/5 {
    generate byte y`t' = rbinomial(1, invlogit(u))
}
quietly reshape long y, i(s) j(t)

meglm y i.t || s: , family(binomial) link(logit) nolrtest nolog

*
* Begin here
*
tempname B V invV v
matrix define `B' = e(b)

matrix define `V' = e(V)
matrix define `V' = `V'[1..`=rowsof(`V')-1', 1..`=colsof(`V')-1']
matrix define `invV' = invsym(`V')

scalar define `v' = `B'[1, `=colsof(`B')']
local column_names : colfullnames(e(b))
local constraint : word `=colsof(`B')' of `column_names'
constraint define 1 _b[`constraint'] = `v'

matrix define `B' = `B'[1, 1..`=colsof(`B')-1']

tempname B1 CD cd
quietly generate double cd = .
quietly levelsof s, local(ss)
foreach s of local ss {
    quietly meglm y i.t if s != `s' || s: , family(binomial) link(logit) constraints(1) nolrtest nolog
    matrix define `B1' = e(b)
    matrix define `B1' = `B1'[1, 1..`=colsof(`B1')-1']
    matrix define `CD' = `B' - `B1'
    matrix define `CD' = `CD' * `invV' * `CD''
    scalar define `cd' = `CD'[1, 1]
    quietly replace cd = `cd' if s == `s'
}
graph twoway dropline cd s if t == 1, msize(small) ylabel( , angle(horizontal) nogrid)

exit

It's in the manner of SAS's PROC MIXED, where it leaves out a cluster at a time, and does not re-fit the random effects covariance matrix.

I think that there are approximations to Cook's distance for generalized linear models that don't require refitting the model, and you might be able to exploit them for mixed effects generalized linear models, too.

(This particular example has a short cut, because the response and predictors are all binary, and so there are only a small number of unique response-predictor patterns that need to be re-fit.)

[Amar Singh]

It works. Thank you very much Joseph.