Simulating new data based on a fitted model object??

[Mark Taper]

Hi,  

I have figured out how to use simulateData specifying the model with a string written with lavaan model syntax.  Now I need to simulate data specifying the model based on a previously fitted model object.  Reading the documentation, I thought that I could do this as simulateData(model=lavaanify(model=MyFittedModel), model.type="sem",sample.nobs=1000).  This gives the error:

Error in as.character.default(x) : 

  no method for coercing this S4 class to a vector

> 

The problem seems to be with lavaanify because   PT=lavaanify(model=MyFittedModel)          gives the error:

Error in as.character.default(x) : 

  no method for coercing this S4 class to a vector

> 

Now I am confused, because MyFittedModel really is a fitted lavaan object.

 summary(MyFittedModel)

lavaan (0.5-17) converged normally after  17 iterations

  Number of observations                          1000

  Estimator                                         ML

  Minimum Function Test Statistic                0.694

  Degrees of freedom                                 4

  P-value (Chi-square)                           0.952

Parameter estimates:

  Information                                 Expected

  Standard Errors                             Standard

                   Estimate  Std.err  Z-value  P(>|z|)

Regressions:

  R ~

    H                 0.406    0.031   12.947    0.000

    LC                0.301    0.031    9.844    0.000

    PC                0.273    0.030    9.042    0.000

    LP                0.225    0.039    5.771    0.000

    SA               -0.063    0.034   -1.832    0.067

    FS               -0.013    0.035   -0.379    0.705

  H ~

    LP                0.386    0.037   10.321    0.000

    LC                0.003    0.031    0.101    0.919

    SA               -0.013    0.030   -0.434    0.664

  LC ~

    LP                0.515    0.035   14.785    0.000

    SA               -0.027    0.035   -0.765    0.444

    FS                0.037    0.033    1.124    0.261

  PC ~

    FS               -0.480    0.029  -16.660    0.000

    LC               -0.002    0.029   -0.060    0.952

  FS ~

    SA                0.517    0.030   17.160    0.000

    LP                0.003    0.034    0.093    0.926

  SA ~

    LP               -0.306    0.034   -9.006    0.000

Variances:

    R                 0.973    0.044

    H                 0.990    0.044

    LC                1.044    0.047

    PC                1.069    0.048

    FS                0.977    0.044

    SA                1.074    0.048

I can't just edit the values into a lavaan syntax string because I need to be able to do this itteratively many times?

I would love any help.

Thanks Mark Taper

[Terrence Jorgensen]

The problem seems to be with lavaanify because   PT=lavaanify(model=MyFittedModel)          gives the error:

Error in as.character.default(x) : 

  no method for coercing this S4 class to a vector

> 

Now I am confused, because MyFittedModel really is a fitted lavaan object.

Right, and the "model" argument of lavaanify() is model syntax (character string) or parameter table (as a list or data.frame).  You can find this in the fitted object:

MyFittedModel@ParTable

It would be nice if you could take the point estimates from the parameterEstimates() output and plug them into the ParTable, then feed that to simulateData()

ptab <- as.data.frame(MyFittedModel@ParTable)

ptab$ustart <- parameterEstimates(MyFittedModel)$est

simulateData(ptab, sample.nobs = 100L)

But that doesn't seem to work.  Perhaps in a future version?...

I can't just edit the values into a lavaan syntax string because I need to be able to do this itteratively many times?

You can, but it is probably way easier to just extract the model-implied moments and simulate from those.

Sigma <- MyFittedModel@Fit@Sigma.hat[[1]]

Mu <- MyFittedModel@Fit@Mu.hat[[1]]

set.seed(123) 

simuData <- rockchalk::mvrnorm(n = 100L, mu = Mu, Sigma = Sigma)

There is a random multinormal generator in MASS, but you can't get the same first 100 observations if you use the same seed and then request N = 101 observations.  The package "rockchalk" revises the function to have this capability, which you may find useful if you do a lot of simulating.

Terry

[Mark L. Taper]

Terry, 

Thanks for your suggestions.  In particular thanks for the idea you say doesn't work:

ptab <- as.data.frame(MyFittedModel@ParTable)

ptab$ustart <- parameterEstimates(MyFittedModel)$est

simulateData(ptab, sample.nobs = 100L)

I looked at ptab, and it still lists most of the parameters as free paramaters.  If you fix them, then the simulation does run.

ptab <- as.data.frame(MyFittedModel@ParTable)

ptab$ustart <- parameterEstimates(MyFittedModel)$est

ptab$free=as.integer(rep(0,dim(ptab)[[1]])) #new code fixing all parameters (I think)

simulateData(ptab, sample.nobs = 100L)

​This solves my problem.  Thanks again for your help.

Best MLT​

-- 
Mark L. Taper:              Environmental & Ecological Analysis

(406) 451-9542

Department of Ecology                 Department of Biology
310 Lewis Hall                               220 Bartram Hall
Montana State University               University of Florida
Bozeman, Montana                        Gainesville, Florida  

Errors like straws upon the surface flow:

Who would search for pearls must dive below.

     -Jon Dryden (1631-1700)

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/yN3u8-Wsirs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at http://groups.google.com/group/lavaan.
For more options, visit https://groups.google.com/d/optout.

[Yves Rosseel]

On 07/15/2015 07:01 PM, Terrence Jorgensen wrote:
> Right, and the "model" argument of lavaanify() is model syntax
> (character string) or parameter table (as a list or data.frame). You
> can find this in the fitted object:
>
> MyFittedModel@ParTable

A safer way is to use the parTable() function to extract the parameter
table, as in

myData <- simulataData( parTable(MyFittedModel) )

> It would be nice if you could take the point estimates from the
> parameterEstimates() output and plug them into the ParTable, then feed
> that to simulateData()
>
> ptab <- as.data.frame(MyFittedModel@ParTable)
> ptab$ustart <- parameterEstimates(MyFittedModel)$est
> simulateData(ptab, sample.nobs = 100L)

Ha! You just got bitten by the nasty as.data.frame() function. The
problem is that -by default- as.data.frame() will convert all strings to
factors, and this confuses the function lavMatrixRepresentation() (at
least in 0.5-18).

Compare these two:

lavMatrixRepresentation(partable = MyFittedModel@ParTable)

lavMatrixRepresentation(partable = as.data.frame(MyFittedModel@ParTable))

So this will work:

ptab <- as.data.frame(MyFittedModel@ParTable, stringsAsFactors = FALSE)

Anyway, I have now included some checks that will convert the 'factor'
columns back to 'character' in lavaan 0.5-19.866

In future versions, the parameter table will also include the 'est'
column, and -if found- this one will be used by simulateData().

Yves.

[Terrence Jorgensen]

So this will work:

ptab <- as.data.frame(MyFittedModel@ParTable, stringsAsFactors = FALSE)

Awesome!  That's really convenient -- thanks Yves!

Terry