Simulation with xtmixed – Stops running

[Matteo Zullo]

Hello everyone,

Someone might be able to help me. The issue has not been addressed in prior posts.

I am running a simulation of a two-level model. The simulation environment is simulate and the multi-level environment is xtmixed. My target is the usual 1,000 reps but the software never reaches that mark, sometimes stopping at 200, sometimes at 300, other times at 10 iterations. It just depends on the seed. I exclude this to be: a) a problem with the simulated dataset because this works fine with regression; b) a problem with xtmixed because Stata is able to run the two-level model outside of the simulation environment.

The simulated outcome is student achievement (score_hat) which is a linear combination of a level-1 variable (IQ), a level-2 variable (urban), a cross-level interaction (urbIQ), a level-1 error (e_ij), and a level-2 error (u_i). I do not show how I create the variables to save some space, but I am doing it correctly.

Thank you.

### Creating simulation
set seed 123
program define montecarlo, rclass
drop _all
tempvar IQ urban schoolID studID u_i e_ij score_hat

# Creating dataset [NOT SHOWN...]

# Running HLM model
quietly xtmixed score_hat IQ urban urbIQ || schoolID:, var reml
return scalar b0HLM=_b[_cons]
return scalar bIQHLM=_b[IQ]
return scalar burbanHLM=_b[urban]
return scalar burbIQHLM=_b[urbIQ]

end


### Requiring 1,000 reps
simulate b0HLM=r(b0HLM) bIQHLM=r(bIQHLM) burbanHLM=r(burbanHLM) burbIQHLM=r(burbIQHLM), reps(1000): montecarlo

[Clyde Schechter]

Well, when you do simulations like this, sometimes a simulated data set will turn out not to be usable. For example, some of your simulated data sets might turn out to have intra-school correlations that are very close to zero, or even negative. Such data sets often cause convergence problems for -mixed- (the correct name for the program formerly known as -xtmixed-, if you are using a recent version of Stata). But I'm just speculating here. The way to find out what is going on is to add the -noisily- option to your -simulate- command. Then you will be able to see for yourself where things are falling apart.

If the problem is what I have guessed, non-convergence of -mixed-, then the way to enable the simulations to proceed is to specify the -iterate()- option in -mixed-, giving some reasonable number of iterations that suffices for most convergent estimations. Then check e(converged) and discard the results from any non-convergent estimation.

[Matteo Zullo]

I specified noisily and I do not get any notification. Stata freezes as before, just with the noisily screen.
On the intra-class correlation: I control it by specifying the distribution of the errors and set it to be 10%. This is a realistic one. I tried with 20% and 30% and the issue persists.

[Joseph Coveney]

Clyde has it. Nevertheless, with such a simple model, mixed should converge with nearly any realistic dataset. I wouldn't have expected to be reliably encountering convergence problems with such a model and any plausible dataset. Perhaps you're not creating the variables correctly? Your declaration of the variables as temporary and yet implementation of them as permanent seems strange, too. Creation of an example dataset takes only eight lines of code in the example down below. I'm not sure what phenomenon it is that you're trying to look into with your simulation, but is your dataset creation so much different that eliding the code saves that much space?

Anyway, you can fit the same model using generalized least squares using

xtset schoolID
xtreg score_hat i.urban##c.IQ, re

which would not give you convergence problems, even in pathological datasets. It doesn't provide restricted maximum likelihood estimates, but you seem to be interested only in the regression coefficient values, anyway, and so it might not matter much for your purposes.

The following goes through the 1000 replications without a hitch. Maybe you can build your datasets off something analogous.

version 15.1

clear *

set seed `=strreverse("1496534")'

program define simem, rclass
    version 15.1
    syntax

    drop _all
    set obs 26
    generate double u = rnormal()
    generate byte sid = _n
    generate byte urb = mod(_n, 2)
    
    expand 100
    generate double iqs = runiformint(80, 120)

    generate double sco = u + (urb - 0.5) / 0.5 + (2 * urb - 1) * (iqs / 100 - 1) + rnormal()

    mixed sco i.urb##c.iqs || sid: , reml dfmethod(kroger) nolrtest nolog

    return scalar urb = _b[1.urb]
    return scalar iqs = _b[iqs]
    return scalar iac = _b[1.urb#c.iqs]
end

simulate iac = r(iac), reps(1000): simem

exit

[Clyde Schechter]

Sorry, I forgot to say that you also need to take the -quietly- off of your -mixed- command so that you can see just what -mixed- is actually doing. At the moment, what is probably the most crucial part is being silenced. I still think the most likely thing is that Stata has not actually frozen but is just cranking iterations on a non-convergent -mixed- estimation. The default number of iterations is 16,000, which can take a very long time to run. When there is no output being shown, it can easily look as if Stata has frozen.

Probably with intra-class correlation at 10, 20 or 30% there won't be the problem I suspected with 0 or negative ICC in a sample. But that isn't the only reason that -mixed- can fail to converge.

Anyway, try it with the noisily option on simulate and the quietly prefix removed from -mixed- and see what happens.

[Joseph Coveney]

With such a simple model, even a deliberately negative ICC doesn't per se faze mixed: run the code below and you'll see that it doesn't balk. Parameterizations of xtreg , re and xtreg , mle allow them to throw out the higher level and set its variance component to zero. In contrast, mixed's default parameterization forces the variance component to remain strictly positive (although, as expected, it collapses); nevertheless, its algorithm reaches convergence by the third iteration.

version 15.1

clear *

set seed `=strreverse("1496557")'

tempname Corr
matrix define `Corr' = J(7, 7, -0.15)
local varlist
forvalues i = 1/`=rowsof(`Corr')' {
    matrix define `Corr'[`i', `i'] = 1
    local varlist `varlist' sco`i'
}

quietly drawnorm `varlist', double corr(`Corr') n(26)
generate byte sid = _n
generate byte urb = mod(_n, 2)

quietly reshape long sco, i(sid) j(pid)
generate int iqs = runiform(80, 120)

xtset sid
xtreg sco i.urb##c.iqs, re
xtreg sco i.urb##c.iqs, mle nolog

mixed sco i.urb##c.iqs || sid: , residuals(exchangeable) noconstant reml nolrtest nolog

// Here
mixed sco i.urb##c.iqs || sid: , reml nolrtest // nolog

exit

[Clyde Schechter]

Joseph Coveney

Thanks, that's interesting. Can you describe what makes a model "simple" for these purposes? I have certainly had the experience of watching -mixed- fail to converge, with the trace showing that one of the lnsig components is marching down to negative infinity (and, of course, never getting there.) Is it just a matter of 2 levels vs more than 2? Or the absence of random slopes? Or what?

[Joseph Coveney]

Clyde, yes, things like that.

There is nothing fancy about it. No quadratic terms, no random slopes, no involved correlation structures fitted to residuals or to random effects. No cross-classification.

The model has only two predictors, one each at two levels, and their interaction. One predictor is categorical with only two categories (and based on what I perceive is the nature of the observational study, neither category is likely to be sparse) and the other predictor is continuous (less likely for near collinearity between the two predictors).

I have certainly had the experience with even what I would consider uncomplicated models to fail to converge, but I wouldn't expect it to be quite so inevitable as the OP sees: every attempt doesn't make it to the end of the simulation run. I can only guess that there is something in the way the datasets are being generated that is liable to give rise to some pathological condition.

[Matteo Zullo]

Thank you all for your reply. I integrated your suggestions in the code posted next. It doesn't run the simulations, but it does run the individual model. I don't know.

version 15.1
clear *
set seed `=strreverse("1496534")'

program define sim, rclass
version 15.1
drop _all
syntax

set obs 20
gen schID = _n
gen double u_i = rnormal(0,10)
gen byte urb = runiform()<0.70

expand 30
bysort schID: gen studID = _n
gen double iq = rnormal(100,15)
gen double iqurb = urb*iq
gen double e_ij = rnormal(0,90)

gen double score = 500 + 0.5*iq + 20*urb + 0.1*iqurb + u_i + e_ij

mixed score iq iqurb i.urb  || schID:, var reml nolrtest

return scalar k =_b[_cons]
return scalar iq = _b[iq]
return scalar urb = _b[1.urb]
return scalar iqurb = _b[iqurb]
end

simulate iqurb = r(iqurb), reps(1000): sim

[Maarten Buis]

It doesn't run the simulations, but it does run the individual model.

If I copy the code you gave us into a .do file and run the entire .do file in one go, then it does run the simulation.

[Matteo Zullo]

Yes, sure, still it stops somewhere before hitting 200.

[Matteo Zullo]

Clyde Schechter Following his handy suggestions, and bracketing for time being the more far-fetched approaches proposed here, I successfully run the 1,000 simulnations with this code. Hope this helps folks in the future. Thank you all.

version 15.1
clear *
set seed `=strreverse("1496534")'
program define sim, rclass
version 15.1
drop _all
syntax

set obs 20
gen schID = _n
gen double u_i = rnormal(0,10)
gen byte urb = runiform()<0.70

expand 30
bysort schID: gen studID = _n
gen double iq = rnormal(100,15)
gen double iqurb = urb*iq
gen double e_ij = rnormal(0,90)
gen double score = 500 + 0.5*iq + 20*urb + 0.1*iqurb + u_i + e_ij

mixed score iq iqurb i.urb  || schID:, var reml nolrtest iterate(1000)
if e(converged)==1{
   return scalar k =_b[_cons]
   return scalar iq = _b[iq]
   return scalar urb = _b[1.urb]
   return scalar iqurb = _b[iqurb]
}
else if e(converged)==0{
estimates clear
}
end
simulate iqurb = r(iqurb), reps(1000): sim

[Joseph Coveney]

bracketing for time being the more far-fetched approaches proposed here

?

I think that most of your problem with occasional failure to converge is that your sample size for the higher level (20 schools) is too small for a residual variance that is targeted to be many times your school variance component. See the summary statistics of the variance component estimates below: the median residual-to-interschool variance ratio is100.

In the extreme, due to the luck of the draw for a substantial proportion of the generated datasets (the average ratio is 3 × 10²²), -mixed- is having trouble teasing them apart with such a small sample size at the higher level of the hierarchy—in one instance where the residual could not be estimated (missing valued), -mixed , reml- put everything in the school random effect.

You can also see the converged results using -mixed , ml- of those four datasets for which -mixed , reml- failed to converge. The estimate for the variance component is zero and that for the residual variance is the better part of 10 000.

In the simulation, you specify a variance for schools of 100 and for students within schools of 8100. If the ratio of the variances is truly of order 100, then you'll need a larger sample size to help assure convergence in those occasional random draws where the ratio happens to be several orders of magnitude. Either that, or use a regularizing prior distribution on the model parameters with a -bayes:- prefix.

.ÿ
.ÿversionÿ15.1

.ÿclearÿ*

.ÿsetÿseedÿ`=strreverse("1496534")'

.ÿ
.ÿglobalÿiterationÿ0

.ÿ
.ÿprogramÿdefineÿsim,ÿrclass
ÿÿ1.ÿÿÿÿÿÿÿÿÿversionÿ15.1
ÿÿ2.ÿÿÿÿÿÿÿÿÿdropÿ_all
ÿÿ3.ÿÿÿÿÿÿÿÿÿsyntax
ÿÿ4.ÿ
.ÿÿÿÿÿÿÿÿÿsetÿobsÿ20
ÿÿ5.ÿÿÿÿÿÿÿÿÿgenÿschIDÿ=ÿ_n
ÿÿ6.ÿÿÿÿÿÿÿÿÿgenÿdoubleÿu_iÿ=ÿrnormal(0,10)
ÿÿ7.ÿÿÿÿÿÿÿÿÿgenÿbyteÿurbÿ=ÿruniform()<0.70
ÿÿ8.ÿ
.ÿÿÿÿÿÿÿÿÿexpandÿ30
ÿÿ9.ÿÿÿÿÿÿÿÿÿbysortÿschID:ÿgenÿstudIDÿ=ÿ_n
ÿ10.ÿÿÿÿÿÿÿÿÿgenÿdoubleÿiqÿ=ÿrnormal(100,15)
ÿ11.ÿÿÿÿÿÿÿÿÿgenÿdoubleÿiqurbÿ=ÿurb*iq
ÿ12.ÿÿÿÿÿÿÿÿÿgenÿdoubleÿe_ijÿ=ÿrnormal(0,90)
ÿ13.ÿ
.ÿÿÿÿÿÿÿÿÿgenÿdoubleÿscoreÿ=ÿ500ÿ+ÿ0.5*iqÿ+ÿ20*urbÿ+ÿ0.1*iqurbÿ+ÿu_iÿ+ÿe_ij
ÿ14.ÿ
.ÿÿÿÿÿÿÿÿÿglobalÿiterationÿ=ÿ$iterationÿ+ÿ1
ÿ15.ÿÿÿÿÿÿÿÿÿmixedÿscoreÿiqÿiqurbÿi.urbÿÿ||ÿschID:,ÿvarÿremlÿnolrtestÿiterate(25)
ÿ16.ÿÿÿÿÿÿÿÿÿifÿ!e(converged)ÿ{
ÿ17.ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿsaveÿ$iteration,ÿreplace
ÿ18.ÿÿÿÿÿÿÿÿÿ}
ÿ19.ÿ
.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿkÿ=_b[_cons]
ÿ20.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿiqÿ=ÿ_b[iq]
ÿ21.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿurbÿ=ÿ_b[1.urb]
ÿ22.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿiqurbÿ=ÿ_b[iqurb]
ÿ23.ÿ
.ÿÿÿÿÿÿÿÿÿtempnameÿB
ÿ24.ÿÿÿÿÿÿÿÿÿmatrixÿdefineÿ`B'ÿ=ÿe(b)
ÿ25.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿsid_uÿ=ÿexp(`B'[1,ÿ6])^2
ÿ26.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿresÿ=ÿexp(`B'[1,ÿ7])^2
ÿ27.ÿÿÿÿÿÿÿÿÿreturnÿscalarÿconvergedÿ=ÿe(converged)
ÿ28.ÿend

.ÿ
.ÿsimulateÿiqurbÿ=ÿr(iqurb)ÿ///
>ÿÿÿÿÿÿÿÿÿsid_uÿ=ÿr(sid_u)ÿresÿ=ÿr(res)ÿconvergedÿ=ÿr(converged),ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿreps(1000):ÿsim

ÿÿÿÿÿÿcommand:ÿÿsim
ÿÿÿÿÿÿÿÿiqurb:ÿÿr(iqurb)
ÿÿÿÿÿÿÿÿsid_u:ÿÿr(sid_u)
ÿÿÿÿÿÿÿÿÿÿres:ÿÿr(res)
ÿÿÿÿconverged:ÿÿr(converged)

Simulationsÿ(1000)
----+---ÿ1ÿ---+---ÿ2ÿ---+---ÿ3ÿ---+---ÿ4ÿ---+---ÿ5ÿ
..................................................ÿÿÿÿ50
..................................................ÿÿÿ100
..................................................ÿÿÿ150
..................................................ÿÿÿ200
..................................................ÿÿÿ250
..................................................ÿÿÿ300
..................................................ÿÿÿ350
..............x...................................ÿÿÿ400
..................................................ÿÿÿ450
..................................................ÿÿÿ500
..................................................ÿÿÿ550
..................................................ÿÿÿ600
..................................................ÿÿÿ650
..................................................ÿÿÿ700
..................................................ÿÿÿ750
..................................................ÿÿÿ800
..................................................ÿÿÿ850
..................................................ÿÿÿ900
..................................................ÿÿÿ950
..................................................ÿÿ1000

.ÿ
.ÿgenerateÿdoubleÿratioÿ=ÿresÿ/ÿsid_u
(1ÿmissingÿvalueÿgenerated)

.ÿsummarizeÿsid_uÿresÿratioÿifÿconverged

ÿÿÿÿVariableÿ|ÿÿÿÿÿÿÿÿObsÿÿÿÿÿÿÿÿMeanÿÿÿÿStd.ÿDev.ÿÿÿÿÿÿÿMinÿÿÿÿÿÿÿÿMax
-------------+---------------------------------------------------------
ÿÿÿÿÿÿÿsid_uÿ|ÿÿÿÿÿÿÿÿ996ÿÿÿÿ115.5572ÿÿÿÿ292.1533ÿÿÿ2.74e-22ÿÿÿÿ8703.03
ÿÿÿÿÿÿÿÿÿresÿ|ÿÿÿÿÿÿÿÿ995ÿÿÿÿÿ8081.17ÿÿÿÿ477.6891ÿÿÿÿ6702.44ÿÿÿ9821.925
ÿÿÿÿÿÿÿratioÿ|ÿÿÿÿÿÿÿÿ995ÿÿÿÿ3.29e+22ÿÿÿÿ8.96e+23ÿÿÿ15.02262ÿÿÿ2.82e+25

.ÿcentileÿsid_uÿresÿratioÿifÿconverged,ÿcentile(25ÿ50ÿ75)

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ--ÿBinom.ÿInterp.ÿ--
ÿÿÿÿVariableÿ|ÿÿÿÿÿÿÿObsÿÿPercentileÿÿÿÿCentileÿÿÿÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+-------------------------------------------------------------
ÿÿÿÿÿÿÿsid_uÿ|ÿÿÿÿÿÿÿ996ÿÿÿÿÿÿÿÿÿ25ÿÿÿÿ13.67959ÿÿÿÿÿÿÿÿ7.071499ÿÿÿÿ21.28644
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ50ÿÿÿÿ82.11065ÿÿÿÿÿÿÿÿ72.76543ÿÿÿÿÿ91.7973
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ75ÿÿÿÿ165.2929ÿÿÿÿÿÿÿÿ154.9986ÿÿÿÿ183.2658
ÿÿÿÿÿÿÿÿÿresÿ|ÿÿÿÿÿÿÿ995ÿÿÿÿÿÿÿÿÿ25ÿÿÿÿ7743.809ÿÿÿÿÿÿÿÿ7710.239ÿÿÿÿ7784.012
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ50ÿÿÿÿ8055.232ÿÿÿÿÿÿÿÿ8027.283ÿÿÿÿ8098.053
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ75ÿÿÿÿ8399.926ÿÿÿÿÿÿÿÿ8351.886ÿÿÿÿ8458.281
ÿÿÿÿÿÿÿratioÿ|ÿÿÿÿÿÿÿ995ÿÿÿÿÿÿÿÿÿ25ÿÿÿÿ48.33989ÿÿÿÿÿÿÿÿ44.62102ÿÿÿÿÿ52.9034
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ50ÿÿÿÿ98.58597ÿÿÿÿÿÿÿÿ86.71426ÿÿÿÿ107.9238
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ75ÿÿÿÿ587.8052ÿÿÿÿÿÿÿÿ376.5247ÿÿÿÿ1174.408

.ÿlistÿsid_u-convergedÿifÿmi(res),ÿnoobsÿabbreviate(20)

ÿÿ+---------------------------+
ÿÿ|ÿÿÿsid_uÿÿÿresÿÿÿconvergedÿ|
ÿÿ|---------------------------|
ÿÿ|ÿ8703.03ÿÿÿÿÿ.ÿÿÿÿÿÿÿÿÿÿÿ1ÿ|
ÿÿ+---------------------------+

.ÿ
.ÿlocalÿfailsÿ:ÿdirÿ"`c(pwd)'"ÿfilesÿ"*.dta"

.ÿ
.ÿforeachÿfailÿofÿlocalÿfailsÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿquietlyÿuseÿ`fail',ÿclear
ÿÿ3.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ_newline(2)ÿ"Iterationÿ`:ÿsubinstrÿlocalÿfailÿ".dta"ÿ""'"
ÿÿ4.ÿÿÿÿÿÿÿÿÿmixedÿscoreÿiqÿiqurbÿi.urbÿÿ||ÿschID:ÿ,ÿnofetableÿnoheaderÿnolog
ÿÿ5.ÿ}


Iterationÿ191

------------------------------------------------------------------------------
ÿÿRandom-effectsÿParametersÿÿ|ÿÿÿEstimateÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
schID:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ2.24e-11ÿÿÿ4.01e-08ÿÿÿÿÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿÿ.
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(Residual)ÿ|ÿÿÿ7776.198ÿÿÿÿ448.961ÿÿÿÿÿÿ6944.212ÿÿÿÿ8707.865
------------------------------------------------------------------------------
LRÿtestÿvs.ÿlinearÿmodel:ÿchibar2(01)ÿ=ÿ1.8e-12ÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ1.0000


Iterationÿ358

------------------------------------------------------------------------------
ÿÿRandom-effectsÿParametersÿÿ|ÿÿÿEstimateÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
schID:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ1.09e-11ÿÿÿ1.16e-10ÿÿÿÿÿÿ8.91e-21ÿÿÿÿ.0133155
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(Residual)ÿ|ÿÿÿ7913.459ÿÿÿ456.8838ÿÿÿÿÿÿÿ7066.79ÿÿÿÿ8861.566
------------------------------------------------------------------------------
LRÿtestÿvs.ÿlinearÿmodel:ÿchibar2(01)ÿ=ÿ0.00ÿÿÿÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ1.0000


Iterationÿ807

------------------------------------------------------------------------------
ÿÿRandom-effectsÿParametersÿÿ|ÿÿÿEstimateÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
schID:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ1.65e-11ÿÿÿ3.66e-08ÿÿÿÿÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿÿ.
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(Residual)ÿ|ÿÿÿ8043.226ÿÿÿ464.3773ÿÿÿÿÿÿ7182.671ÿÿÿÿ9006.884
------------------------------------------------------------------------------
LRÿtestÿvs.ÿlinearÿmodel:ÿchibar2(01)ÿ=ÿ0.00ÿÿÿÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ1.0000


Iterationÿ969

------------------------------------------------------------------------------
ÿÿRandom-effectsÿParametersÿÿ|ÿÿÿEstimateÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
schID:ÿIdentityÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(_cons)ÿ|ÿÿÿ5.68e-11ÿÿÿ7.13e-08ÿÿÿÿÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿÿ.
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿvar(Residual)ÿ|ÿÿÿ8366.164ÿÿÿ483.0264ÿÿÿÿÿÿÿ7471.05ÿÿÿÿ9368.523
------------------------------------------------------------------------------
LRÿtestÿvs.ÿlinearÿmodel:ÿchibar2(01)ÿ=ÿ0.00ÿÿÿÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ1.0000

.ÿ
.ÿexit

endÿofÿdo-file


.

[Matteo Zullo]

Thank you for your thorough explanation. The final paper will use, as the case of most recent literature, real estimates coming from a real dataset.

[Joseph Coveney]

Thank you for your thorough explanation. The final paper will use, as the case of most recent literature, real estimates coming from a real dataset.

Well, okay, but if your real estimates coming from a real dataset have a disparity in the variance components as great as what you show above, then you'll benefit from one of those far-fetched approaches proposed in this thread, namely, to reparameterize the mixed model from one of variance components to its exactly equivalent one of compound symmetric repeated measures. The model parameter estimates are identical; the reparameterized variance components are also exactly equivalent (and easily interconverted afterward if you're interested in them).

I've run your do-file below exactly as you wrote it, with just the model reparameterization and a convergence flag to prove its benefit—100% convergence for the thousand repetitions. See below (the changes indicated by "<- here".)

.ÿ
.ÿversionÿ15.1

.ÿclearÿ*

.ÿsetÿseedÿ`=strreverse("1496534")'

.ÿ
.ÿprogramÿdefineÿsim,ÿrclass
ÿÿ1.ÿversionÿ15.1
ÿÿ2.ÿdropÿ_all
ÿÿ3.ÿsyntax
ÿÿ4.ÿ
.ÿsetÿobsÿ20
ÿÿ5.ÿgenÿschIDÿ=ÿ_n
ÿÿ6.ÿgenÿdoubleÿu_iÿ=ÿrnormal(0,10)
ÿÿ7.ÿgenÿbyteÿurbÿ=ÿruniform()<0.70
ÿÿ8.ÿ
.ÿexpandÿ30
ÿÿ9.ÿbysortÿschID:ÿgenÿstudIDÿ=ÿ_n
ÿ10.ÿgenÿdoubleÿiqÿ=ÿrnormal(100,15)
ÿ11.ÿgenÿdoubleÿiqurbÿ=ÿurb*iq
ÿ12.ÿgenÿdoubleÿe_ijÿ=ÿrnormal(0,90)
ÿ13.ÿ
.ÿgenÿdoubleÿscoreÿ=ÿ500ÿ+ÿ0.5*iqÿ+ÿ20*urbÿ+ÿ0.1*iqurbÿ+ÿu_iÿ+ÿe_ij
ÿ14.ÿ
.ÿ*ÿmixedÿscoreÿiqÿiqurbÿi.urbÿÿ||ÿschID:,ÿvarÿremlÿnolrtest
.ÿmixedÿscoreÿiqÿiqurbÿi.urbÿÿ||ÿschID:,ÿresiduals(exchangeable)ÿnoconstantÿremlÿnolrtestÿ//ÿ<-ÿhere
ÿ15.ÿ
.ÿreturnÿscalarÿkÿ=_b[_cons]
ÿ16.ÿreturnÿscalarÿiqÿ=ÿ_b[iq]
ÿ17.ÿreturnÿscalarÿurbÿ=ÿ_b[1.urb]
ÿ18.ÿreturnÿscalarÿiqurbÿ=ÿ_b[iqurb]
ÿ19.ÿreturnÿscalarÿconvergedÿ=ÿe(converged)ÿ//ÿ<-ÿhere
ÿ20.ÿend

.ÿ
.ÿsimulateÿiqurbÿ=ÿr(iqurb)ÿ///
>ÿÿÿÿÿÿÿÿÿconvergedÿ=ÿe(converged)ÿ///ÿ<-ÿhere
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ,ÿreps(1000):ÿsim

ÿÿÿÿÿÿcommand:ÿÿsim
ÿÿÿÿÿÿÿÿiqurb:ÿÿr(iqurb)
ÿÿÿÿconverged:ÿÿe(converged)

Simulationsÿ(1000)
----+---ÿ1ÿ---+---ÿ2ÿ---+---ÿ3ÿ---+---ÿ4ÿ---+---ÿ5ÿ
..................................................ÿÿÿÿ50
..................................................ÿÿÿ100
..................................................ÿÿÿ150
..................................................ÿÿÿ200
..................................................ÿÿÿ250
..................................................ÿÿÿ300
..................................................ÿÿÿ350
..................................................ÿÿÿ400
..................................................ÿÿÿ450
..................................................ÿÿÿ500
..................................................ÿÿÿ550
..................................................ÿÿÿ600
..................................................ÿÿÿ650
..................................................ÿÿÿ700
..................................................ÿÿÿ750
..................................................ÿÿÿ800
..................................................ÿÿÿ850
..................................................ÿÿÿ900
..................................................ÿÿÿ950
..................................................ÿÿ1000

.ÿ
.ÿtabulateÿconverged

e(convergedÿ|
ÿÿÿÿÿÿÿÿÿÿ)ÿ|ÿÿÿÿÿÿFreq.ÿÿÿÿÿPercentÿÿÿÿÿÿÿÿCum.
------------+-----------------------------------
ÿÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿÿ1,000ÿÿÿÿÿÿ100.00ÿÿÿÿÿÿ100.00
------------+-----------------------------------
ÿÿÿÿÿÿTotalÿ|ÿÿÿÿÿÿ1,000ÿÿÿÿÿÿ100.00

.ÿ
.ÿexit

endÿofÿdo-file


.