You are here

Errors when implementing FML function

4 posts / 0 new
Last post
Jannik's picture
Offline
Joined: 10/01/2018 - 04:42
Errors when implementing FML function
AttachmentSize
Binary Data FML.R2.43 KB

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

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
why?

OpenMx does the transformation automatically in the C++ code. Why do you want to do it yourself?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
some advice
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.

Jannik's picture
Offline
Joined: 10/01/2018 - 04:42
Thank you!

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.