Attachment | Size |
---|---|
FML.R [6] | 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