Error in Sixway(), something with RL?

Welcome to the forum for R2MLwiN users. Feel free to post your question about R2MLwiN here. The Centre for Multilevel Modelling take no responsibility for the accuracy of these posts, we are unable to monitor them closely. Do go ahead and post your question and thank you in advance if you find the time to post any answers!

Go to R2MLwiN: Running MLwiN from within R >> http://www.bris.ac.uk/cmm/software/r2mlwin/
Post Reply
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Error in Sixway(), something with RL?

Post by GerineLodder »

I get the following error when I try to run the sixway function:
Error in if (RL1$resmatrix[1] == "Error") { :
missing value where TRUE/FALSE needed
In addition: Warning message:
In par(new = TRUE) : calling par(new=TRUE) with no plot
It doesn't occur when I run the analyses with 100 iterations, but when I use 1000 or more the problem does occur.

Perhaps good to know: I run the sixway function from another function (because I have many models and it saves time):
( for: model, I fill in the model name)

Code: Select all

printSixway <- function(model, name = "model", width = 500, height = 500, ...) {
  modName <- paste("sixway_","_",sep = name)
  modName2 <- paste("trajectories_", "", sep = name)
  modName3 <- paste("","_", sep = name)
  
  trajectories(model)
  #dev.print(png,paste(modName2,".png"), width = width, height = height)
  
    for (i in colnames(model["chains"])) {
    sixway(model["chains"][,i],i, name = paste(modName3, "", sep=i))
    #dev.print(png,paste(modName,".png",sep = i),
    #width = width, height = height)
    }
}
ChrisCharlton
Posts: 1351
Joined: Mon Oct 19, 2009 10:34 am

Re: Error in Sixway(), something with RL?

Post by ChrisCharlton »

To help diagnose this could you please run the following on the problematic model and report the output?:

Code: Select all

# Load coda library for raftery.diag function
library(coda)
# Combine multiple chains into one
flatchain <- as.matrix(model@chains)
# Manually run Raftery-Lewis diagnostic
RL1 <- raftery.diag(flatchain, q = 0.025, r = 0.005, s = 0.95, converge.eps = 0.001)
# Examine results
RL1$resmatrix
str(RL1)
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Re: Error in Sixway(), something with RL?

Post by GerineLodder »

I did that, and it seems the problem may lie in the variance of the intercept at level 1?
> RL1$resmatrix
M N Nmin I
deviance 25 30105 3746 8.04
FP_Intercept 12 12718 3746 3.40
FP_predictor 18 19852 3746 5.30
RP3_var_Intercept 126 156618 3746 41.80
RP2_var_Intercept 120 111108 3746 29.70
RP1_var_bcons_1 NA NA 3746 NA
> str(RL1)
List of 2
$ params : Named num [1:3] 0.005 0.95 0.025
..- attr(*, "names")= chr [1:3] "r" "s" "q"
$ resmatrix: num [1:6, 1:4] 25 12 18 126 120 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:6] "deviance" "FP_Intercept" "FP_predictor" "RP3_var_Intercept" ...
.. ..$ : chr [1:4] "M" "N" "Nmin" "I"
- attr(*, "class")= chr "raftery.diag"
These are the model results:
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.01) multilevel model (Binomial)
N min mean max N_complete min_complete mean_complete max_complete
mlwinschool 99 7 43.77778 210 99 7 42.98990 206
mlwinclass 234 1 18.52137 36 232 1 18.34483 36
Estimation algorithm: MCMC Elapsed time : 45.99s
Number of obs: 4256 (from total 4334) Number of iter.: 10000 Chains: 1 Burn-in: 10000
Bayesian Deviance Information Criterion (DIC)
Dbar D(thetabar) pD DIC
4359.628 4293.372 66.255 4425.883
---------------------------------------------------------------------------------------------------
The model formula:
logit(vicGLdC3, denom) ~ 1 + predictor + (1 | mlwinschool) + (1 |
mlwinclass)
Level 3: mlwinschool Level 2: mlwinclass Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Cred. Interval] ESS
Intercept -1.07761 0.08627 -12.49 8.3e-36 *** -1.24652 -0.91409 952
predictor -0.29973 0.10732 -2.79 0.005226 ** -0.50867 -0.08293 834
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------------------------------------------
The random part estimates at the mlwinschool level:
Coef. Std. Err. [95% Cred. Interval] ESS
var_Intercept 0.07757 0.03478 0.01881 0.15397 166
---------------------------------------------------------------------------------------------------
The random part estimates at the mlwinclass level:
Coef. Std. Err. [95% Cred. Interval] ESS
var_Intercept 0.07061 0.05936 0.00698 0.20360 12
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err. [95% Cred. Interval] ESS
var_bcons_1 1.00000 1e-05 1.00000 1.00000 10000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
And this is the model specification:

Code: Select all

flog01_vicGLdC3 = logit(vicGLdC3,denom) ~ 1 + predictor + (1|mlwinschool) + (1|mlwinclass)
ChrisCharlton
Posts: 1351
Joined: Mon Oct 19, 2009 10:34 am

Re: Error in Sixway(), something with RL?

Post by ChrisCharlton »

Yes, it does appear to be the bcons parameter that is causing the problem as the NA returned from the raftery.diag() function will give NA in the comparison that the message quoted (RL1$resmatrix[1] == "Error"):

Code: Select all

> NA == "Error"
[1] NA
So the check should be:

Code: Select all

if (is.na(RL1$resmatrix[1]))
Running the raftery.diag() function directly gives the same message:

Code: Select all

> raftery.diag(mymodel1@chains[, "RP1_var_bcons_1"])

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
Error in if (x$resmatrix[1] == "Error") cat("\nYou need a sample size of at least",  : 
  missing value where TRUE/FALSE needed
so I wonder if something has changed in R since this was written. In any case you should be able to work around it by checking your "i" value to see if it is "RP1_var_bcons_1", and if so skipping the call to sixway() for this parameter.
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Re: Error in Sixway(), something with RL?

Post by GerineLodder »

Hi Chris,

Could you perhaps indicate how to do that?
Isn't there a problem in the output as well, given that the estimates for the constant are also odd?

Gerine
ChrisCharlton
Posts: 1351
Joined: Mon Oct 19, 2009 10:34 am

Re: Error in Sixway(), something with RL?

Post by ChrisCharlton »

The following change should do this:

Code: Select all

printSixway <- function(model, name = "model", width = 500, height = 500, ...) {
  modName <- paste("sixway_","_",sep = name)
  modName2 <- paste("trajectories_", "", sep = name)
  modName3 <- paste("","_", sep = name)
  trajectories(model)
  #dev.print(png,paste(modName2,".png"), width = width, height = height)
  for (i in colnames(model["chains"])) {
    if (i != "RP1_var_bcons_1") {
      sixway(model["chains"][,i],i, name = paste(modName3, "", sep=i))
      #dev.print(png,paste(modName,".png",sep = i), width = width, height = height)
    }
  }
}
As the bcons.1 chain will always only contain a series of the value one then the MCMC diagnostics don't make sense for this parameter.
GerineLodder
Posts: 42
Joined: Thu Jul 28, 2016 11:33 am

Re: Error in Sixway(), something with RL?

Post by GerineLodder »

thanks, worked like a charm!
Post Reply