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"):

So the check should be:

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!