MCMC starting values

[NilsGYork]

Hello

(- sorry, not sure in which forum this fits best -)

I'm trying to estimate a multilevel mixed response model in MLwiN using MCMC and own starting values (all run through runmlwin). Unfortunately, MLwiN crashes every time. I use the example from the runmlwin help file to illustrate my problem ('Two-level mixed bivariate continuous and binary response probit model (fitted using IGLS MQL1)'). I only get a Windows error message, not one from MLwiN.

Code: Select all

* Load data

use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear
generate binlrt = (standlrt>0)

* Run IGLS model (works fine!)	

runmlwin 	(normexam cons, equation(1)) ///
(binlrt cons, equation(2)) ///
, level1(student:(cons, equation(1)) (cons, equation(2))) ///
		discrete(distribution(normal binomial) link(probit) denominator(cons cons)) ///
		nosort nopause batch mlwinpath("C:\Program Files (x86)\MLwiN v2.28\i386\mlwin.exe") maxiterations(250)

* Run MCMC with IGLS starting values (works fine!)	

runmlwin 	(normexam cons, equation(1)) ///
		(binlrt cons, equation(2)) ///
		, level1(student:(cons, equation(1)) (cons, equation(2))) ///
		discrete(distribution(normal binomial) link(probit) denominator(cons cons)) ///
		nosort nopause batch  mlwinpath("C:\Program Files (x86)\MLwiN v2.28\i386\mlwin.exe") ///
		mcmc(on burnin(50) chain(200)) initsprevious

* Specify own starting values
// Betas		
matrix BB = [0,0,1,0.5,1]
			
// Covariance matrix with last row set to zeros
	matrix A = J(2,2,0.125)
	matrix B = J(3,3,0.125)
	matrix C = J(2,3,0)
	matrix A = A, C
	matrix B = C', B
	matrix V = A \ B
	forvalues i = 1/5 {
		matrix V[`i',`i'] = 0.125
		matrix V[5,`i'] = 0
	}
mata: st_replacematrix("V", makesymmetric(st_matrix("V")))			

* Run MCMC with own starting values (doesn't work)			
runmlwin 	(normexam cons, equation(1)) ///
		(binlrt cons, equation(2)) ///
		, level1(student:(cons, equation(1)) (cons, equation(2))) ///
		discrete(distribution(normal binomial) link(probit) denominator(cons cons)) ///
		nosort nopause batch  mlwinpath("C:\Program Files (x86)\MLwiN v2.28\i386\mlwin.exe") ///
		mcmc(on burnin(50) chain(200)) initsb(BB) initsv(V)	

I tried this with a bivariate linear model and it worked fine. So I assume this has something to do with constraining the variance of the probit error term to 1?

The matrices BB and V look like this:

Code: Select all

BB[1,5]
    c1  c2  c3  c4  c5
r1   0   0   1  .5   1

Code: Select all

symmetric V[5,5]
      c1    c2    c1    c2    c3
r1  .125
r2  .125  .125
c1     0     0  .125
c2     0     0  .125  .125
c3     0     0     0     0     0

Cheers
Nils

[ChrisCharlton]

Thanks for letting us know about this. We will look into the cause and hope to fix this in a later release of the software. Incidentally, I found that if I used the same structure of starting values as you but changed the values of 0.125 into 0.00125 it did run successfully, so the crash may relate to your choice of values.

[NilsGYork]

Hello Chris,
thanks for your reply.

Well spotted! I hadn't noticed this behaviour before because I had only played around with values larger than 0.125. I now checked and it seems that the problem seems to be with the variance of the covariance term (V[4,4]). If this is chosen to be very small, the other terms on the diagonal of the variance-covariance matrix can be bigger. But I guess specifying a small variance for the covariance term as starting value would imply highly informative priors, right? Not really what I want.

So I'm wondering if this is a problem with MLwiN or my choice of starting values...

[billb]

Hi Nils,
Have been discussing your emails with Chris and we will look at why the starting values cause a crash and get back to you. It may be that your starting values are generating an initial position for the variance function which implies a negative variance for certain observations which would cause the routine not to be able to be started - as I say I'll look into this.
HOWEVER in these models it is important not to get confused between starting values and priors. For models with complex level 1 variances we essentially use a uniform prior over the parameter space of the parameters in the variance function and so there is no link between starting values used and prior distribution. Hope this makes sense?
Regards,
Bill.