Login

News

Errors when implementing FML function

Posted on
No user picture. Jannik Joined: 10/01/2018
Forums
Fit Functions
Hello,
I’m trying to implement a slightly changed Maximum Likelihood function using mxFitfunctionR and ran into some problems that I wasn’t able to solve. Here is a short version of my current attempt at implementing a covariance-based FML function without any changes:

simple_FML <- function(model, state) {
obsCov <- model$data$observed# observed covariance
nmanif <- ncol(obsCov) # number of observed variables

Sigma_try_out <- try(mxGetExpected(model,"covariance")) # expected covariance matrix

if (class(Sigma_try_out)== "try-error"){return(NA)} # if impossible, return NA
# if it can be built, I am calculating the FML:
else {
inv_Sigma_try <- try(solve(Sigma_try_out))
if (class(inv_Sigma_try)== "try-error"){return(NA)}
else{
FML_temp <- log(det(Sigma_try_out))+tr(obsCov%*%inv_Sigma_try)-log(det(obsCov))- nmanif
if (anyNA(model$data$means)){return(FML_temp)} #this FML_temp is reported if there are no means in the model
else {
obsMeans <- model$data$means
M_try <- try(mxGetExpected(model,"means")) # expected means
if (class(M_try)== "try-error" ){return(NA)}
else{FML_temp <- FML_temp + (obsMeans- M_try)%*%inv_Sigma_try%*% t(obsMeans- M_try) #compute FML with means
return(FML_temp)}}
}
}
}

The template for this FML function can be found here: http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_calis_sect074.htm

The resulting parameter values are quite close to the real ones – still they are constantly off. And the standard errors are nowhere near the ones reported when using the default fitfunction. Unfortunately I have no idea why this is the case.
I attatched a file where I simulated some data and then used the default fitfunction and the one reported above. Does anyone know what went wrong?

If there is a nice workaroud for obtaining the FIML and FML results for a current model with an OpenMx function, that would be even better than my approach. I only want to add a term to the normal FML or FIML value. I already tried the following, but it didn't work:

basefitfunction <- mxFitFunctionML()

my_fitfunction <- function(model, state){
fit <- mxEval(model = model, expression = basefitfunction, compute = TRUE)
return(fit)
}

my_fitfunction <- mxFitFunctionR(my_fitfunction)

myModel <- mxModel(baseModel, my_fitfunction)
mxRun(myModel)

I would very much appreciate your help.
Jannik

Replied on Thu, 10/04/2018 - 11:19
Picture of user. jpritikin Joined: 05/23/2012

OpenMx does the transformation automatically in the C++ code. Why do you want to do it yourself?
Log in or register to post comments
Replied on Thu, 10/04/2018 - 12:16
Picture of user. AdminRobK Joined: 01/24/2014

The resulting parameter values are quite close to the real ones – still they are constantly off. And the standard errors are nowhere near the ones reported when using the default fitfunction.

I tried running your script in this environment,


> mxVersion()
OpenMx version: 2.11.3 [GIT v2.11.3-dirty]
R version: R version 3.4.0 (2017-04-21)
Platform: i386-w64-mingw32
Default optimiser: NPSOL
NPSOL-enabled?: Yes
OpenMP-enabled?: No

, and the largest absolute discrepancy between point estimates is on the order of 1e-3. I'd say they're pretty darn close. The standard errors from using your fitfunction are too large by a factor of sqrt(1000), because the Hessian at the solution needs to be scaled by sample size to get the resulting standard errors, which OpenMx does automatically when using an MxFitFunctionML.

Note that your -2logL at the solution will necessarily differ from OpenMx's, because OpenMx uses a normalized loglikelihood (i.e., OpenMx does not drop constants from the loglikelihood).

Your fitfunction is also not production-ready because it has no contingency to deal with invertible but indefinite model-expected covariance matrices. See this other thread for details. Instead of inverting with solve(), try first Choleky-factoring the expected covariance matrix with chol(), and then solving for the inverse from the Cholesky factor with chol2inv(). Be sure to wrap chol() in try() and have your function return `NA` in the case of an error. Inverting via Cholesky factorization is also more numerically accurate, which might reduce the discrepancies between the point estimates from using your fitfunction versus the built-in fitfunction.

I already tried the following, but it didn't work:

Unfortunately, fitfunctions cannot be evaluated via mxEval(). Fitfunctions have to be evaluated via a call to the backend (e.g., mxRun()). But, you don't want to call the backend from an R fitfunction. Per the documentation for mxFitFunctionR(), "fitfun should not call R functions that use OpenMx's compiled backend, including (but not limited to) `omxMnor()`, because doing so can crash R."

I only want to add a term to the normal FML or FIML value

If that's all you want to do, then try putting your MxModel object named "Model" inside another MxModel object. Then, use an MxFitFunctionAlgebra to define the container MxModel's fitfunction to be `Model.fitfunction` plus whatever extra term you want to add to it.

Log in or register to post comments
Replied on Fri, 10/05/2018 - 08:09
No user picture. Jannik Joined: 10/01/2018

Thank you very much for your help! Using the "nested" Models, as you suggested, is a far better way than my current approach - especially because it is faster than my own FIML implementation.
Log in or register to post comments
  Login

News