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

FML.R | 2.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

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

I tried running your script in this environment,

, 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.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."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.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.