You are here

Estimating a non-linear model in OpenMx

6 posts / 0 new
Last post
pgok's picture
Offline
Joined: 11/13/2019 - 16:40
Estimating a non-linear model in OpenMx

Hi, I'm pretty new to OpenMx, having only used it off-and-on over the past few years. That said, I am pretty interested in using it more, especially given the multilevel modeling capabilities (e.g., as discussed in Pritkin et al., 2017). However, my question isn't about the multilevel modeling capabilities, my question is about fitting non-linear models.

I am looking at trying to fit an exponential decay/growth model. It would look something like this: y~b2+(b0-b2)*exp(b1*time). Where b1 is constrained to be negative (so that everything after b2 eventually drops out of the model as time goes to infinity). The problem I'm running into is that it isn't clear to me how to include a parameter in the exponential part of the model. I ran across this thread on including non-linear constraints, but I can't seem to figure out how to adapt that code to my issue. A (not working) minimal example is below:

  time=seq(0, 100, by=0.01)
  y<-1+(5-1)*exp(-.1*time)+rnorm(length(time),sd=.2)
 
  data_list<-data.frame(
    time=time,
    y=y
  )
 
  plot(time,y)
 
  modelExp<-mxModel("exp", type="RAM", manifestVars=c(names(data_list)),
                    latentVars = "b1",
                    mxPath(from="b1", to="b1", free=F, value=0, arrows=2),
                    mxPath(from="one", to="b1", value=-.1),
                    mxAlgebra(exp(b1*data.time), name="expPart"),
                    mxPath(from="one", to="y", label="b2"),
                    mxPath(from="time", to="y", label="expPart"),
                    mxPath(from="one", to="time"),
                    mxData(data_list, type="raw")
                    )
  fit<-mxRun(modelExp)

Part of the problem is that I can't figure out how to reference the parameter b1 (which I am using a zero variance latent variable for) in the mxAlgebra. The other part of the problem is, based on the previous thread I mentioned, I don't see a way to estimate the (b0-b2) part. I don't think I actually need to estimate those as two separate parameters (I already estimate b2 as a kind of intercept, so I can just estimate b0-b2 as some single parameter b* and then do the algebra separately), but it looks like the model with try to make "expPart" equal to b0-b2.

It's entirely possible that I am either trying to do something entirely unreasonable, or something trivially easy in an overly complicated way, but regardless it escapes me how I would do this and any help would be greatly appreciated.

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Nonlinear Regression

Hi!

Here's how I would do the model you proposed.

  modelExp<-mxModel("exp", type="RAM", manifestVars='y',
                    mxMatrix('Full', nrow=1, ncol=3,
                      labels=c('b0', 'b1', 'b2'), name='params',
                      free=TRUE, values=c(.5, -1, 1)),
                    mxAlgebra(b2+(b0-b2)*exp(b1*data.time), name='predY'),
                    mxPath(from="one", to="y", labels='predY[1,1]', free=FALSE),
                    mxPath(from="y", arrows=2, labels='resid', values=1),
                    mxData(data_list, type="raw")
                    )
  fit<-mxRun(modelExp)

It seems to do okay.

 
> summary(fit)
Summary of exp 
 
free parameters:
   name matrix row col    Estimate    Std.Error A
1 resid      S   y   y  0.04087680 0.0005780564  
2    b0 params   1   1  5.00507642 0.0130146711  
3    b1 params   1   2 -0.09956287 0.0005180620  
4    b2 params   1   3  0.99726107 0.0026133824
pgok's picture
Offline
Joined: 11/13/2019 - 16:40
Thanks that does seem to work

Thanks that does seem to work! Huge help.

pgok's picture
Offline
Joined: 11/13/2019 - 16:40
Follow up question: is it

Follow up question: is it possible to include a latent variable in the mxAlgebra expression? For example, if I had some latent variable "L" that I wanted to include in an arbitrary expression [say, y=b1*log(b2*L)], I am guessing that there is a way to do that, but it doesn't seem to have a clear analog to the manifest variable case [i.e., there isn't a way to do b1*log(b2*latent.L)].

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Not exactly

A latent variable doesn’t have a value per se. it is a variable with a distribution. One possibility would be to create a mixture distribution of say 12 latent trait values weighted by Gaussian quadrature weights. This is currently one way of handling such things, but it doesn’t scale well as with 2 latent moderators it would take 144 evaluations to get the likelihood of a data vector, and generally with m quadrature points and d latent moderators, m^d gets large fast. There may be an alternative in OpenMx 3.0, but it’s not ready yet. Some time in 2022 looks likely though. The other route would be via MCMC, but there appear to be many pitfalls and limitations in that area.

pgok's picture
Offline
Joined: 11/13/2019 - 16:40
Thanks for the response. I

Thanks for the response. I can't say I'm terribly surprised that it gets radically more complex for the latent variable case, but I figured I'd ask anyways. I'll be interested to see what OpenMx 3.0 has in store for this though. Thanks again!