Multi level mixed effect model with Inverse Gaussian family and log link function

[ashar malik]

Hello, I am trying to run a two level model with Inverse Gaussian family and log link function using meglm. But it does not support this specification. My specification in Stata is as follows
meglm depmvar indepvar || clusvar: , family(inversegaussian) link(log)

next I tried gsem command but same problem. Installed runmlwin but it don't support inverse Gaussian.
I used Gaussian invgaussian but none working. IS there any alternative command or user written program that i can use to fix the problem.

[Joseph Coveney]

. . . I am trying to run a two level model with Inverse Gaussian family and log link function using meglm. But it does not support this specification. . . . IS there any alternative command or user written program that i can use to fix the problem.

You could use bayesmh with a custom log-likelihood and flat or weak priors. I show how below. Begin at the "Begin here" comment; the stuff above is to create a toy dataset for illustration.

version 19

clear *

// seedem
set seed 688788178

// Upper level (clusters)
quietly set obs 50
generate byte pid = _n
generate double pid_u = rnormal()

// Lower level (individual observations)
quietly expand 10
generate double pre = runiform(-0.5, 0.5)
generate double xbu = 1 + 2 * pre + pid_u
generate double out = rigaussian(exp(xbu), 1)

*
* Begin here
*
glm out c.pre, family(igaussian) link(log) vce(cluster pid) nolog

#delimit ;
bayesmh out , 
    likelihood(llf(-((out - exp({out: pre U[pid], xb}))^2 / out / 
        exp({out: pre U[pid], xb})^2 + 3 * ln(out) + ln(2 * _pi)) / 2)) 
    prior({out:}, flat) prior({var_U}, halfcauchy(0, 10)) 
    initial({out:pre} 0 {out:_cons} 1 {var_U} 1) 
    nomodelsummary nodots;
#delimit cr

exit

A few notes:

1. I first also fit a marginal (population average) model using glm for comparison. It didn't recover the regression parameters (1.0 for the intercept and 2.0 for the slope of the continuous predictor) so accurately as the hierarchical / multilevel model fitted with bayesmh, which got the coefficients spot-on.

2a. For production work, you'll probably want to run at least two chains with dispersed initial values. And if you want to beautify the diagnostic graphical evaluations (bayesgraph) you can play around with longer chains or, equivalently, thinning.

2b. Although I didn't do it above, I recommend putting regularizing priors on at least the regression coefficients, even if you choose to keep a flat prior on the intercept.

3. You might be able to smarten up the nonlinear substitutable expression for the custom log-likelihood by use of the define() option, if esthetics are that important.

4. An alternative is to code a maximum-likelihood regression model in Mata using an instance of its Quadrature() class to integrate over the random effects. I've done that for a linear mixed model just as an exercise, and using its default adaptive Gauss-Kronrod integration method gave results identical to the corresponding meglm command with its default adaptive quadrature method. My exercise fitted a two-level model as you describe (i.e., a single random intercept), and it converged lickety-split.

5. If Stata doesn't implement a particular combination of distribution family and link function, then there's usually a reason. Although the example above went without a hitch, I'm guessing that this multilevel generalized linear model might be finicky to fit to less-than-pretty real datasets—good luck!

Complete do-file and log file attached.

[ashar malik]

Thank you, very valuable. The only difference is that the "pre" variable in my case is binary: requesting to specify the level for the latent variable "pre'"

[ashar malik]

While running your above syntax, it gives "distribution halfcauchy is not supported in opion prior() .

[Joseph Coveney]

The only difference is that the "pre" variable in my case is binary: requesting to specify the level for the latent variable "pre'"

Well, if your binary variable is coded 0/1, then you don't have to do anything—just include it as-is.

I'm not sure why you are referring to it as a "latent variable". Predictors in these kinds of model are manifest variables and not latent factors.

While running your above syntax, it gives "distribution halfcauchy is not supported in opion prior() .

What version of Stata are you using?

[ashar malik]

Stata 18, sorry, there was an error in the message. The correct message is "distribution halfcauchy is not supported in option prior()".

[ashar malik]

Please review the log. I have appended the log file you shared.

[ashar malik]

When I run this for my model, the glm runs, but for bayesmh it gives the error "no level specified for latent variable pre"

[Joseph Coveney]

Stata 18, sorry, there was an error in the message. The correct message is "distribution halfcauchy is not supported in option prior()".

See here. You'll need either to upgrade or to choose an alternative distribution. Depending upon how much data (how many clusters) you have, you might be able to get away with a flat prior on this parameter estimate, too.

When I run this for my model, the glm runs, but for bayesmh it gives the error "no level specified for latent variable pre"

See here.