Wednesday, 29 Mar 2006 Lars Korsholm <[email protected]> asked about how to
use -xtmixed- to fit the model with heterogeneous errors w_j:
> I'll like to fit a model like
>
> X_ij = a_i + e_ij, i = 1...20, j = 1,2,3
>
> where a_i are iid normal whith mean = a and sd = sigma and where e_ij are iid
> normal with mean =0 and sd=w_j, i.e. different residual error on different
> days
>
> I have tried with xtmixed, but I cant make it work, any hint?
There is no direct way to use -xtmixed- to fit multilevel models with
heterogeneous errors. Lars can try -gllamm- to solve his problem. For example,
. tabulate day, gen(daycat)
. eq het: daycat1 daycat2 daycat3
. gllamm x, i(id) s(het)
where day defines days (j=1,2,3), id defines subjects (i=1,...,20) and x
defines the dependent variable.
It is worth mentioning, however, that in some cases one can use -xtmixed- to
fit models with heterogeneous errors. For example, if for the model above we
assume that the heterogeneity is of the form, for example, Var(eij) =
w(j)*sigma_e^2, i.e. residual variance is proportional to the constant we can
use -xtmixed- to estimate unknown parameters sigma and sigma_e. w(j) is some
covariate that varies across the levels of day and residual variance is
proportional to its values.
The model may be expressed as follows:
X_ij = a + u_i + w_j*u_ij, i = 1...20, j = 1,2,3 (2)
where a is the overall mean, w_j are known for each day =1,2,3, u_i's are iid
N(0,sigma^2), and u_ij's are iid N(0, sigma_e^2). Disregarding w_j, model (2)
is the simplest one-level model as defined in [XT] xtmixed with u_i's defining
the first level and u_ij's defining errors. Alternatively we can treat u_ij's
as another level defining an interaction between day and subject (subject =
1,...,20). Then we have a two-level model with random intercept model for the
first level, slope-only model for the second level and zero variability for
the residual error.
Let's now look at the example. The data in the example is simulated according
the model (2) with a = 5, sigma = 2 and sigma_e = 0.08.
/***************************** begin ****************************************/
clear
set seed 12345
local a = 5
local sigma = 2
local sigma_e = 0.08
set obs 3
gen day = _n
gen w = uniform()
expand 20
sort day, stable
by day: gen id = _n
sort id, stable
by id: gen double u = (_n==1)*`sigma'*invnorm(uniform())
by id: replace u = sum(u)
gen double y = `a' + u + w*`sigma_e'*invnorm(uniform())
xtmixed y || id: || day: w, nocons
/****************************** end *****************************************/
The output from -xtmixed- is the following:
------------------------------------------------------------------------------
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: log restricted-likelihood = 25.133721 (not concave)
Iteration 1: log restricted-likelihood = 26.869516 (not concave)
Iteration 2: log restricted-likelihood = 28.568708 (not concave)
Iteration 3: log restricted-likelihood = 28.930586 (not concave)
Iteration 4: log restricted-likelihood = 29.222204
Iteration 5: log restricted-likelihood = 29.797077
Iteration 6: log restricted-likelihood = 29.834883
Iteration 7: log restricted-likelihood = 30.065664
Iteration 8: log restricted-likelihood = 30.08004
Iteration 9: log restricted-likelihood = 30.081975 (not concave)
Iteration 10: log restricted-likelihood = 30.08207
Iteration 11: log restricted-likelihood = 30.082082
Iteration 12: log restricted-likelihood = 30.082082
Computing standard errors:
Mixed-effects REML regression Number of obs = 60
-----------------------------------------------------------
| No. of Observations per Group
Group Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
id | 20 3 3.0 3
day | 60 1 1.0 1
-----------------------------------------------------------
Wald chi2(0) = .
Log restricted-likelihood = 30.082082 Prob > chi2 = .
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 4.523483 .451127 10.03 0.000 3.63929 5.407676
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity |
sd(_cons) | 2.017478 .3272741 1.468005 2.772619
-----------------------------+------------------------------------------------
day: Identity |
sd(w) | .0831976 .0093019 .0668255 .1035809
-----------------------------+------------------------------------------------
sd(Residual) | 4.62e-06 9.35e-06 8.77e-08 .0002435
------------------------------------------------------------------------------
LR test vs. linear regression: chi2(2) = 313.05 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference
------------------------------------------------------------------------------
The 95% CIs for sigma = sd(_cons) (1.468005, 2.772619) and sigma_e = sd(w)
(.0668255, .1035809) contain the true parameter values 2 and 0.08,
respectively.
-- Yulia
[email protected]
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
