discrete multivariate models

[iandohoo]

Hi
Is it possible to fit multivariate models to discrete data using -runmlwin- (eg. multivariate logistic models)?
Thank you
Ian Dohoo

[GeorgeLeckie]

Hi Ian,

Yes you can fit multivariate multilevel multivariate probit model (but not logit model)...

Code: Select all

* Open the tutorial data
use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear

* Generate a unique student identifier
drop student
gen student = _n

* Dichotmize normexam
generate binexam = (normexam>0)

* Dichotmize standlrt
generate binlrt = (standlrt>0)

* Fit two-level bivariate response probit model by quasilikelihood methods
runmlwin ///
	(binexam cons, equation(1)) ///
	(binlrt cons, equation(2)) ///
	, ///
	level2(school: (cons, eq(1)) (cons, eq(2))) ///
	level1(student: (cons, eq(1)) (cons, eq(2))) ///
	discrete(distribution(binomial binomial) link(probit) denominator(cons cons)) nopause

* Fit two-level bivariate response probit model by MCMC methods
runmlwin ///
	(binexam cons, equation(1)) ///
	(binlrt cons, equation(2)) ///
	, ///
	level2(school: (cons, eq(1)) (cons, eq(2))) ///
	level1(student: (cons, eq(1)) (cons, eq(2))) ///
	discrete(distribution(binomial binomial) link(probit) denominator(cons cons)) ///
	mcmc(on) initsprevious nopause

Best wishes


George

[iandohoo]

Hi George
Thank you for the code for fitting multivariate probit models. Everything is working fine for me until I try to constrain one of the coefficients to be constant over the responses.
Here is the code and error message. The outcome is repeated measurements of "highscc" (binary) and there are 2 predictors "c_heifer" and "t_dim".

. runmlwin (highscc0 cons c_heifer t_dim0, eq(1)) ///
> (highscc1 cons c_heifer t_dim1, eq(2)) ///
> (highscc2 cons c_heifer t_dim2, eq(3)) ///
> (c_heifer, eq(1/3)), ///
> level2(herdid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3))) ///
> level1(cowid: (cons, eq(1)) (cons, eq(2)) (cons, eq(3))) ///
> discrete(distribution(binomial binomial binomial) ///
> denominator(cons cons cons) ///
> link(probit)) nopause

Incorrect number of distributions specified.

The code works fine if I omit the line (c_heifer, eq(1/3)) ... but I think that is what I need to add to get a constant coefficient for that predictor.

Am I missing something obvious.

Thanks
Ian Dohoo

[GeorgeLeckie]

Hi Ian,

Thanks for you post. You have found a bug. We were checking the number of distributions/denominators against the number of fixed part brackets rather than the number of model equations. We have now implemented a fix and this will be available in the next released of runmlwin which will be available in a few weeks. We will email you the fix off-forum so that you can carry on with your modelling.

Below we provide an example where we constrain the coefficient of girl to be constant across the two response equations of a two-level (students-within-schools) bivariate binary response (binexam and binlrt) probit model.

Best wishes

George


The commands are

Code: Select all

* Load the data
use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear

* Generate a unique student identifier
drop student
gen student = _n

* Dichotmize normexam
generate binexam = (normexam>0)

* Dichotmize standlrt
generate binlrt = (standlrt>0)

* Fit two-level bivariate response probit model by quasilikelihood methods
runmlwin ///
   (binexam cons, equation(1)) ///
   (binlrt cons, equation(2)) ///
   (girl, equation(1/2)) ///
   , ///
   level2(school: (cons, eq(1)) (cons, eq(2))) ///
   level1(student: (cons, eq(1)) (cons, eq(2))) ///
   discrete(distribution(binomial binomial) link(probit) denominator(cons cons)) nopause

The commands with output are

Code: Select all

. * Load the data
. use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear

. 
. * Generate a unique student identifier
. drop student

. gen student = _n

. 
. * Dichotmize normexam
. generate binexam = (normexam>0)

. 
. * Dichotmize standlrt
. generate binlrt = (standlrt>0)

. 
. * Fit two-level bivariate response probit model by quasilikelihood methods
. runmlwin ///
>    (binexam cons, equation(1)) ///
>    (binlrt cons, equation(2)) ///
>    (girl, equation(1/2)) ///
>    , ///
>    level2(school: (cons, eq(1)) (cons, eq(2))) ///
>    level1(student: (cons, eq(1)) (cons, eq(2))) ///
>    discrete(distribution(binomial binomial) link(probit) denominator(cons cons)) nopause
 
MLwiN 2.27 multilevel model                     Number of obs      =      4059
Multivariate response model
Estimation algorithm: IGLS, MQL1

-----------------------------------------------------------
                |   No. of       Observations per Group
 Level Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
         school |       65          2       62.4        198
-----------------------------------------------------------

Run time (seconds)   =       4.34
Number of iterations =          5
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
binexam      |
      cons_1 |  -.0670705   .0595326    -1.13   0.260    -.1837523    .0496113
-------------+----------------------------------------------------------------
binlrt       |
      cons_2 |  -.0787224   .0503415    -1.56   0.118    -.1773898    .0199451
-------------+----------------------------------------------------------------
     girl_12 |   .1584247    .044779     3.54   0.000     .0706595    .2461899
------------------------------------------------------------------------------

------------------------------------------------------------------------------
   Random-effects Parameters |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
Level 2: school              |
                 var(cons_1) |   .1563027   .0327351      .0921431    .2204623
          cov(cons_1,cons_2) |   .0811783   .0219956      .0380676     .124289
                 var(cons_2) |   .0918002   .0211727      .0503026    .1332979
-----------------------------+------------------------------------------------
Level 1: student             |
                var(bcons_1) |          1   7.89e-12             1           1
        cov(bcons_1,bcons_2) |   .4227481   .0119733      .3992808    .4462154
                var(bcons_2) |          1   1.37e-10             1           1
------------------------------------------------------------------------------

[Paatrick]

Hi George

Has that problem really been fixed in the current version? Because I think I ran into the same problem:

Code: Select all

runmlwin (total_cholesterol       cons, equation(1)) ///
         (systolic_blood_pressure cons, equation(2)) ///
         (age, equation(1/2)), ///
         level1(id: (cons, equation(1)) (cons, equation(2))) ///
         nosort

gives me "Incorrect number of distributions specified." My runmlwin has the Distribution-Date: 20130324 and Iam using Stata 12.

Thanks for your consideration!
Patrick

Edit/P.S.: I should mention that my outcomes are both continuous, not discrete.

[GeorgeLeckie]

Dear Patrick,

Apologies we have still not put the latest version of runmlwin on SSC

I will endeavor to do this soon

Meanwhile I will email you the fix off forum

Best wishes

George