Problems with convergence when parameters are fixed and question about BIC

6 replies [Last post]
metavid's picture
Offline
Joined: 06/23/2014

Good day,

I am having some trouble with posting due to the spam filter being continuously triggered so I also provided an attachment for the R code in addition to the data set.

I have two different questions.
My questions revolve around multiple group factor analysis (3 groups). There are a lot more groups but for the sake of the questions and computation speed, I provided the code for three groups.
Thank you for your suggestions and comments.

1. Fixing of model parameters: I fixed all the parameters in a one factor model. The code works for one group (model1) but when three groups are put into one entire run of code (mgroups), I get the warning: "Error: The job for model 'new' exited abnormally with the error message: Objective function returned a value of NaN at iteration 0.1." How could I circumvent this issue to let the code work for the three groups in one run?

2. Total of BIC from one entire code run unequal with the sum of BIC from individual models
I would like to compute the information criteria using three groups in one whole code. I specified the code for three groups in one run (newmodel) and the code for three separate groups (m1,m2,m3). I further specified the code below to compute the information criteria. The thing is, the sum of the AIC's from three separate model is exactly equal with the AIC computed from one entire code run (newmodel) but it is not the case for the BIC's. They do not correspond at all. Is the BIC computed differently when multiple groups are specified in one code or is it just ordinary to find that the sum of the BIC's from separate groups does not equal the BIC when the three groups are combined in one code?

AttachmentSize
codemx.txt9.03 KB
dataset.txt11.71 KB
AdminRobK's picture
Offline
Joined: 01/24/2014
Oddly enough, I also

Oddly enough, I also triggered the spam filter when I tried to reply with my non-admin account.

Concerning your second question: in general, BIC isn't additive. For a given model and dataset, BIC = -2logL + log(N)*p, where p = number of free parameters. If we have Model_1 and Model_2, each applied to separate and independent datasets, and having no free parameters in common,
BIC_1 = -2logL_1 + log(N_1)*p_1
BIC_2 = -2logL_2 + log(N_2)*p_2
If we combine the two models into a supermodel,
overall BIC = -2logL_1 + -2logL_2 + log(N_1 + N_2)*(p_1 + p_2) != BIC_1 + BIC_2

Concerning your first question: I had to make a few small edits to your syntax to reproduce the error you're reporting. Even so, there's nothing obviously wrong with your setup, so I'm not sure why optimization fails on its first function evaluation for you. You could try the following syntax, which uses the rarely-invoked independent=T argument, and see if it does what you want:

mgroups <- mxModel("mg",
                   mxModel("g1",independent = T,
                           mxMatrix( type="Full", nrow=14, ncol=1, values=rep(1,14), free=rep(F,14),name="facLoadingsE82" ),
                           mxMatrix( type="Symm", nrow=1, ncol=1, values=1, free=F, name="facVariancesE82" ),
                           mxMatrix( type="Diag", nrow=14, ncol=14, values=rep(1,14), free=rep(F,14), name="resVariancesE82" ),
                           mxMatrix( type="Full", nrow=1, ncol=14, values=rep(3.5,14), free=rep(F,14), name="varMeansE82" ),
                           mxMatrix( type="Full", nrow=1, ncol=1, values=1, free=F,name="facMeansE82" ),
                           mxAlgebra( expression= facLoadingsE82 %&% facVariancesE82 + resVariancesE82,name="expCovE82" ),
                           mxAlgebra(expression= varMeansE82 + t(facLoadingsE82 %*% facMeansE82),name="expMeanE82" ),
                           mxData(traintE[traintE$year==1982,],type="raw"),
                           mxFIMLObjective(covariance="expCovE82", means="expMeanE82", dimnames=extn)
                   ),
                   mxModel("g2",independent = T,
                           mxMatrix( type="Full", nrow=14, ncol=1, values=rep(1,14), free=rep(F,14),name="facLoadingsE83" ),
                           mxMatrix( type="Symm", nrow=1, ncol=1, values=1, free=F, name="facVariancesE83" ),
                           mxMatrix( type="Diag", nrow=14, ncol=14, values=rep(1,14), free=rep(F,14),name="resVariancesE83" ),
                           mxMatrix( type="Full", nrow=1, ncol=14, values=rep(3.5,14), free=rep(F,14),name="varMeansE83" ),
                           mxMatrix( type="Full", nrow=1, ncol=1, values=1, free=F,name="facMeansE83" ),
                           mxAlgebra( expression= facLoadingsE83 %&% facVariancesE83 + resVariancesE83,name="expCovE83" ),
                           mxAlgebra(expression= varMeansE83 + t(facLoadingsE83 %*% facMeansE83),name="expMeanE83" ),
                           mxData(traintE[traintE$year==1983,],type="raw"),
                           mxFIMLObjective(covariance="expCovE83", means="expMeanE83", dimnames=extn)
                   ),
                   mxModel("g3",independent = T,
                           mxMatrix( type="Full", nrow=14, ncol=1, values=rep(1,14), free=rep(F,14),name="facLoadingsE84" ),
                           mxMatrix( type="Symm", nrow=1, ncol=1, values=1, free=F, name="facVariancesE84" ),
                           mxMatrix( type="Diag", nrow=14, ncol=14, values=rep(1,14), free=rep(F,14), name="resVariancesE84" ),
                           mxMatrix( type="Full", nrow=1, ncol=14, values=rep(3.5,14), free=rep(F,14) ,name="varMeansE84" ),
                           mxMatrix( type="Full", nrow=1, ncol=1, values=1, free=F,name="facMeansE84" ),
                           mxAlgebra( expression= facLoadingsE84 %&% facVariancesE84 + resVariancesE84,name="expCovE84" ),
                           mxAlgebra(expression= varMeansE84 + t(facLoadingsE84 %*% facMeansE84),name="expMeanE84" ),
                           mxData(traintE[traintE$year==1984,],type="raw"),
                           mxFIMLObjective(covariance="expCovE84", means="expMeanE84", dimnames=extn)),
                   mxAlgebra(g1.objective + g2.objective + g3.objective,
                             name="LL"),
                   mxAlgebraObjective("LL")
)
mgroupsrun <- mxRun(mgroups)

EDIT: The multigroup model above will run even without the independent=T argument under the OpenMx 2.0 beta.

metavid's picture
Offline
Joined: 06/23/2014
Summary statistics

Good day,

I was able to run the code above, using independent=T. The thing now is, that I have to have the summary statistics, specifically the information criteria (AIC and BIC with both parameter and df penality). mxRun went well but now I can't seem to compute the information criteria given by summary() due to a warning: argument 'invSDs' is missing, with no default. Is there a way to resolve this or should I stick with the formulae, given that the LL for the model is provided from mxRun (but not the AIC/BIC)?

Thanks.

AdminRobK's picture
Offline
Joined: 01/24/2014
Try summary(mgroupsrun$g1), etc.

Sorry, I should have mentioned that summary() works differently with independent=TRUE included in every submodel. In such a case, each submodel is handled separately from all the others, and the computer doesn't bundle them together as parts of a supermodel. You're really only getting a -2loglikelihood for the supermodel because you requested it as an algebra.

You can see a summary for each submodel with summary(mgroupsrun$g1), and likewise for g2 and g3. Since you know the formulae for BIC with different penalties, thanks to M. Hunter below, you could just use the supermodel -2loglikelihood and compute the supermodel BIC from it.

Using independent=TRUE shouldn't be necessary to do what you want to do in the first place, so the problem that prompted you to make your original post is unexpected, and is probably a bug--fortunately one not present in the 2.0 beta.

metavid's picture
Offline
Joined: 06/23/2014
Another question

Hi,

Thank you for your elaborate response.

I have another question.
1. Are the formulae you provided (for BIC) merely applicable to the BIC from Openmx with parameter penalty?
2. How does the formula for the BIC taking into account degrees-of-freedom penalty look like?

Again, thank you.

mhunter's picture
Offline
Joined: 07/31/2009
Applies to both

1. The formulas were given for the parameters penalty, but the same idea applies for the df penalty. Neither penalty has they additive property you're looking for.

Parameters penalty
BIC_1 = -2logL_1 + log(N_1)*p_1
BIC_2 = -2logL_2 + log(N_2)*p_2

overall BIC = -2logL_1 + -2logL_2 + log(N_1 + N_2)*(p_1 + p_2) != BIC_1 + BIC_2

2. Here's the df penalty.
df penalty
BIC_1 = -2logL_1 + log(N_1)*df_1
BIC_2 = -2logL_2 + log(N_2)*df_2

overall BIC = -2logL_1 + -2logL_2 + log(N_1 + N_2)*(df_1 + df_2) != BIC_1 + BIC_2

The key is that the log(N_1 + N_2) can't be broken apart.

metavid's picture
Offline
Joined: 06/23/2014
Appreciate it

Thanks again for the clarification!