Hi metaSEM users,

As I have posted previously in this forum, I am using metaSEM to conduct a network meta-analysis. I have been attempting to fit a metaSEM model using both ML and REML estimation using some sample data from a meta-analysis of the effectiveness of diabetes interventions. I was able to successfully conduct the analysis using ML estimation, but for some reason when I attempt to obtain a fit using the reml() function, I obtain the following error:

Running Variance component with REML

Error in running the mxModel:

Error in reml(trt, S, RE.constraints = con) :

trying to get slot "runstate" from an object (class "simpleError") that is not an S4 object

In addition: Warning message:

The job for model 'Variance component with REML' exited abnormally with the error message: Objective function returned an infinite value at iteration 0.1.

What's particularly strange is that I have been able to successfully obtain ML and REML model fits using these same data using the mvmeta routine in STATA. Moreover, I am able to exactly replicate the results from the ML fit using the meta() function in metaSEM in STATA. So, I don't think that the error I observed is due to some weirdness in the data.

I have attached the code and data that generated this error. Is it possible that there is a bug in reml() that is responsible for this error?

Thanks in advance for your help.

Kind regards,

- Patrick

Attachment | Size |
---|---|

diabetes analysis for metSEM forum.R | 2.06 KB |

diabetes4.csv | 1.15 KB |

Based on this thread I tried using the value of Tau (the estimated variance of the study effect sizes) I obtained from mvmeta in STATA as a starting value which I passed to RE.startvalues (for reference, the value of Tau from mvmeta in STATA is 0.0007673). Unfortunately, I'm still getting the same error message that I documented above.

- Patrick

Hi Patrick,

I confirm that I also have the same error message. As a matter of fact, this error also occurs in some of my data. The implementations of the REML estimation method in reml() and reml3() in the metaSEM package were based on Cheung (2013). The idea is to transform the data and the model so that the fixed-effects are removed from the parameters before the estimation. There are only one row with lots of columns of variables and lots of constraints. Because of the constraints, numerical stability might be an issue.

I have not figured out how to solve it yet. I hope that others may have ideas on how to improve the numerical stability of this approach.

BTW, you coded the missing values as 0 in your example. Did you intend to do so?

Cheung, M. W.-L. (2013). Implementing restricted maximum likelihood estimation in structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 20(1), 157–167. doi:10.1080/10705511.2013.742404

Mike

Hi Mike,

Thanks so much for replying! Do you have any leads on what might be the cause of this error? For me, it only seems to occur after the mxRun statement in your code for the reml function, but I'm having a hard time figuring out how to continue debugging after that. I also don't receive an error for all the test data that I use (see attached for an example that does not throw an error).

I'll try reading the paper that you suggested in your comments (thanks for providing that reference, by the way).

As for your question about the handling of missing values, bear with me as I give a slightly off-topic explanation. The input data for network meta-analysis is a series of comparisons between a reference treatment and a series of other treatments. For example, say you're doing a network meta-analysis on the effectiveness of various diabetes treatments A, B, C, and D. If you make A the reference treatment, you would need to calculate the treatment comparisons AB, AC, and AD for each study. If you wish to view the other comparisons not directly represented in this scheme (e.g., the comparison BC), you can derive that comparison by taking advantage of the consistency equations, which state that the comparisons not directly modeled in the meta-analysis are linear combinations of the comparisons that are directly modeled in the meta-analysis (for more information, see Salanti, 2012; attached).

However, what happens if some of the values for the reference treatment (in this case, treatment A) for a given study are missing? In this case, it seems like you would be unable to calculate the AB, AC, and AD comparisons for that study, and therefore would be unable to include that study in the meta-analysis, even though that study presumably includes some of the other comparisons that you're interested in (e.g., the BC comparisons).

In this situation, the recommendation of White, Barrett, Jackson, and Higgins (2012; attached) is to represent the unobserved treatment arms with a small amount of information (i.e., if your effect size is the standardized mean difference, represent the unobserved arms with a small mean and a large standard deviation). Taking this approach, the calculated comparisons involving the reference treatment will be relatively unaffected, but you can still use the consistency equations to utilize the information from the focal study about, for example, the BC comparisons. This is the approach that I'm using in the example script that I attached to my first message.

Sorry, I did actually notice that you gave some details about how numerical stability may be part of the issue, but on re-reading the post I realized that it looked like I hadn't read your post. By "do you have any idea what the cause of the error is", I meant, do you know what part of the metaSEM and / or OpenMx code is generating the error? I'd like to help figure out how to tackle this issue, but I'm a little uncertain how to even begin the process of debugging in this example.

The error message comes from OpenMx. It means that after the metaSEM package generated an OpenMx model and started trying to find the best values for its free parameters, OpenMx computed the fit function and found a NA or NAN result. This bad result happened at major iteration 0, minor iteration 1. This typically happens with poor starting values. You could try adjusting the starting values, and/or checkpointing the optimization using mxOptions.

Following Michael's suggestion, one may prefer to run the model manually. This may give a better control on the mxOptions. For example,

mod.reml <- remlNotRun(trt, S, RE.constraints = con)

mod.fit <- mxRun(mod.reml)

Thanks very much for the suggestions, Mike and Michael! I'll try this out tomorrow morning.

Hey Mike, quick question, I see in the function you attached that you've added some code that uses the Diag() function. I don't seem to have any packages with the Diag function. Could you please let me know which package to install to obtain Diag()?

diag, not Diag is in base R.

OpenMx has chosen to not implement this function in the backend.

`?mxAlgebra`

Lists the implemented functions that are legal in mxAlgebras

These include diag2vec and vec2diag

Diag() is available in metaSEM version 0.8-4. diag(c("a","b")) works in R 2.15 but not in R 3.0.0. You may refer to http://r.789695.n4.nabble.com/Behaviors-of-diag-with-character-vector-in-R-3-0-0-td4663735.html for the changes. Since metaSEM needs to create diagonal matrices of vectors, I created a version that works for characters.

I have installed the latest version of metaSEM, which fixed the problem with Diag(). Unfortunately, so far the checkpointing has not been particularly enlightening for me; the code

mod <- remlNotRun(trt, S, RE.constraints = con)

mod.run <- mxRun(mod, checkpoint = TRUE)

results in the following error:

"This application has requested the Runtime to terminate it in an unusual way. Please contact the application's support team for more information."

After which R crashes.

Meanwhile, the .omx file created by R only contains the following information:

iterations timestamp objective a

0 "Oct 01 2013 10:10:52 AM" 1.#INF0 0.10000

Which only suggests to me that whatever error is occurring takes effect immediately in the estimation process.

As I have mentioned further up in the thread, I have actually changed the start values already, including obtaining the correct variance component start values from a different program (Stata) and attempting to run this model with those start values. So far, changing the start values has not helped.