Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
st: Separate intercept mixed model
From
Michael Mitchell <[email protected]>
To
[email protected]
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
(432 real changes made)
. * create 2 groups, trt and control
. generate trt = (id > 24)
. * make weight 10 pounds more for treatment
. replace weight = weight + 10 if trt == 1
(216 real changes made)
. * 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:
Performing gradient-based 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:
Performing gradient-based 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/