Possible bug in R2MLwiN

[autarkie]

Comparing the output for my model in MLwiN and R, I have found serious discrepancies. In R, estimates were printed as belonging to different variables, even though the numbers themselves were identical. Here is a re-run of my model with `debugmode=TRUE`.

The R version is as follows:

Code: Select all

-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN multilevel model (Normal) 
Estimation algorithm:  MCMC      Elapsed time : 829.38s 
Number of obs:  16673 (from total 18841 )        Number of iter.: 20000          Burn-in: 10000 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
115547.648 115506.289 41.353     115589.000 
--------------------------------------------------------------------------------------------------- 
The model formula:
unemp100 ~ (0 | cons + educationgm + HT + gender + StartCareer60 + 
    StartCareer60:HT + avejobgm + avejobgm:HT + carlengthgm) + 
    (1 | cons + StartCareer60 + StartCareer60:HT) + (2 | cons + 
    StartCareer60 + StartCareer60:HT)
Level 2: country     Level 1: MERGEID      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                      Coef.   Std. Err.      z    Pr(>|z|)       [95% Cred.   Interval]     ESS 
cons                3.60620     0.52371   6.84   7.826e-12  ***     2.55728     4.65196   10000 
educationgm        -0.14639     0.01522  -9.61   7.384e-22  ***    -0.17663    -0.11658    9372 
HT                  1.48281     0.21762   6.78    1.24e-11  ***     1.06275     1.91989   10000 
gender              1.01506     0.12152   8.38   5.144e-17  ***     0.78305     1.25600   10000 
StartCareer60      -0.10891     0.02355  -4.61   4.033e-06  ***    -0.15567    -0.06341   10000 
StartCareer60:HT   -1.60088     0.61555  -2.62    0.008886  **     -2.82614    -0.40949    7913 
avejobgm           -0.14369     0.02075  -6.94   3.941e-12  ***    -0.18491    -0.10401   10668 
avejobgm:HT         0.13507     0.04394   3.10    0.001946  **      0.05060     0.22456    9714 
carlengthgm         0.77653     0.24015   3.25    0.001134  **      0.30106     1.24803   10000 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  

And the MLwiN output is this:



Notice how the values are switched. The estimate of -1.601 related to "avejobgm" variable in MLwiN (which is actually correct) is related to StartCareer60:HT interaction in R output. The value of -0.144 pertaining to "carlengthgm" in the MLwiN output corresponds to "avejobgm" in R output. Here is the R call I use:

Code: Select all

mlmod <- runmw(unemp100~(0|cons+educationgm+HT+gender+StartCareer60+StartCareer60:HT + avejobgm + avejobgm:HT + carlengthgm) + (1|cons+StartCareer60+StartCareer60:HT) + (2|cons+StartCareer60+StartCareer60:HT),
               levID=c("country","MERGEID"),data=ml.data,
               estoptions=list(EstM=1,debugmode=T,clre=covmatrix,mcmcOptions=list(hcen=2),
                               mcmcMeth=list(iterations=20000,burnin=10000,thinning=2,Lev1VarM=3)))

I don't do anything else to the model, just call it from R using `debugmode=TRUE` option, take a screenshot in MLwiN, close MLwiN, and then see model output in R by calling it via the object name. There is a warning in MLwiN saying that my variable names may cause problem in Stata (which I don't use). To me it looks like the estimates are sent back to R in the proper order (if you compare R and MLwiN output without looking at variable names, the order of the estimates is identical), but the order of the variable names is changed. The interaction terms are grouped together after main terms in MLwiN, whereas R retains the order from the original formula.

Additional info: MLwiN version is 2.28, R2MLwiN version is 0.1-8, R version is 3.0.2.

[richardparker]

Hi - many thanks for highlighting this. It appears to occur whenever an interaction is specified in the R2MLwiN formula (as variable1:variable2) - R2MLwiN places this at the end of the list of predictors as modelled in MLwiN, regardless of where it was specified in the original formula, but then doesn't take this into account when importing the results back in.

Whilst we devise a solution in R2MLwiN, a workaround in the meantime is to separately create your interactions in R (rather than expressing them as variable1:variable2 in the R2MLwiN formula).

Best wishes,

Richard

[richardparker]

Hi - the latest version of R2MLwiN (0.2-0), released today, now has a bug fix for this.

Thanks again for reporting it.

Best wishes,

Richard