Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

# st: Separate intercept mixed model

 From Michael Mitchell To statalist@hsphsun2.harvard.edu Subject st: Separate intercept mixed model Date Tue, 8 Feb 2011 00:58:41 -0800

```Greetings

I am trying to estimate a mixed (multilevel) model using a separate
intercept approach, and am getting an error saying that the likelihood
evaluates to missing. Here is the model I am estimating with the
error...

. webuse pig, clear
(Longitudinal analysis of pig weights)
. replace week = week - 1
. * create 2 groups, trt and control
. generate trt = (id > 24)
. * make weight 10 pounds more for treatment
. replace weight = weight + 10     if trt == 1
. * Model 0:
. * Enter trt as factor variable
. xtmixed weight ibn.trt c.week , nocons || id:

Performing EM optimization:
likelihood evaluates to missing
r(430);

I try to model this using a conventional coding for the variable
-trt-, and all works well...

. * Model 1:
. * Predict weight from week and treatment
. * express treatment as a factor variable
. xtmixed weight i.trt c.week || id:

Performing EM optimization:

Iteration 0:   log restricted-likelihood = -1015.4632
Iteration 1:   log restricted-likelihood = -1015.4632

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =       432
Group variable: id                              Number of groups   =        48

Obs per group: min =         9
avg =       9.0
max =         9

Wald chi2(2)       =  25363.93
Log restricted-likelihood = -1015.4632          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.trt |         11   1.144154     9.61   0.000     8.757499     13.2425
week |   6.209896   .0390633   158.97   0.000     6.133333    6.286458
_cons |   25.06551     .82399    30.42   0.000     23.45052     26.6805
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
sd(_cons) |    3.90138   .4198202      3.159527    4.817419
-----------------------------+------------------------------------------------
sd(Residual) |   2.096356   .0757444      1.953034    2.250195
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   470.28 Prob >= chibar2 = 0.0000

And I am able to get the fitted mean by treatment group for week 0, shown below.

. * show means by treatment
. margins trt , at(week=0)

Adjusted predictions                              Number of obs   =        432

Expression   : Linear prediction, fixed portion, predict()
at           : week            =           0

------------------------------------------------------------------------------
|            Delta-method
|     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trt |
0  |   25.06551     .82399    30.42   0.000     23.45052     26.6805
1  |   36.06551     .82399    43.77   0.000     34.45052     37.6805
------------------------------------------------------------------------------

Now, I try modeling this using a separate intercept model, but fitting
it using manually constructed dummy variables and it fits nicely. I
feel like this model is identical to model 0, but model 0 fails and
model 2 succeeeds. Further estimates of the parameters -trt0- and
-trt1- from model 2 match the estimates from the previous margins
command from model 1.

. * Model 2:
. * Now express trt as a dummy variable
. * Use a separate intercept model, getting
. * estimates of the mean weight for week
. * 0 by treatment group
. gen trt0 = (trt==0)
. gen trt1 = (trt==1)
. xtmixed weight trt0 trt1 week, nocons || id:

Performing EM optimization:

Iteration 0:   log restricted-likelihood = -1015.4632
Iteration 1:   log restricted-likelihood = -1015.4632

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =       432
Group variable: id                              Number of groups   =        48

Obs per group: min =         9
avg =       9.0
max =         9

Wald chi2(3)       =  34743.66
Log restricted-likelihood = -1015.4632          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trt0 |   25.06551     .82399    30.42   0.000     23.45052     26.6805
trt1 |   36.06551     .82399    43.77   0.000     34.45052     37.6805
week |   6.209896   .0390633   158.97   0.000     6.133333    6.286458
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
sd(_cons) |    3.90138   .4198202      3.159527    4.817419
-----------------------------+------------------------------------------------
sd(Residual) |   2.096356   .0757444      1.953034    2.250195
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   470.28 Prob >= chibar2 = 0.0000

. * the estimate of the difference in means matches trt from model 1
. lincom trt1 - trt0

( 1)  - [weight]trt0 + [weight]trt1 = 0

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
(1) |         11   1.144154     9.61   0.000     8.757499     13.2425
------------------------------------------------------------------------------

So, I cannot figure out how to fit model 2 while expressing -trt-
using factor variable notation, i.e. ibn.trt, rather than manually
fitting dummy variables, i.e. -trt0- and -trt1-.

When I fit model 0, omitting the random effects portion of the
model, it seems that the estimate for treatment group 0 is being
omitted even though I specified -ibn.trt- and -nocons-, see below.

. xtmixed weight ibn.trt c.week , nocons

Mixed-effects REML regression                   Number of obs      =       432

Wald chi2(2)       =         .
Log restricted-likelihood =          .          Prob > chi2        =         .

------------------------------------------------------------------------------
weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trt |
0  |  (omitted)
1  |   22.39341   1.078251    20.77   0.000     20.28008    24.50675
|
week |    9.62792   .1601441    60.12   0.000     9.314043    9.941796
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
sd(Residual) |   12.74731   .4336723      11.92505    13.62627
------------------------------------------------------------------------------

Thank you for your patience in reading this far! I would be most
grateful for suggestions.

Many thanks,

Michael Mitchell
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```