Page 1 of 1

Error in Sixway(), something with RL?

Posted: Tue May 01, 2018 11:09 am
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)
    }
}

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

Posted: Wed May 02, 2018 9:06 am
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)

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

Posted: Fri May 04, 2018 9:49 am
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)

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

Posted: Fri May 04, 2018 10:15 am
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.

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

Posted: Tue May 08, 2018 12:52 pm
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

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

Posted: Tue May 08, 2018 1:02 pm
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.

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

Posted: Thu May 17, 2018 9:49 am
by GerineLodder
thanks, worked like a charm!