simulateData - lavaan Error: wrong number of arguments in modifier

[Ash Gillis]

Hi there,

I've fitted a multi-group mediation model with three parallel mediators and three covariates.

The model successfully converged and provides a summary just fine.

However, when I attempt to simulate the data using simulateData or using sim in the package simsem, I get the following error:

Error in lavaanify(model = model, meanstructure = meanstructure, int.ov.free = int.ov.free,  : 

  lavaan ERROR: wrong number of arguments in modifier (M1YDirectA,M1YDirectB,M1YDirectC,M1YDirectD) of element Y~M1

However, the number of modifiers I've included is equal to the number of groups I have, which is 4. 

I saw a previous post where someone had the same issue and the suggestion was to try it as a multilevel model. Is there a different way to do this or perhaps my construction of the model is incorrect? I constructed the model in two different ways because I am unsure about whether each indirect and total effect needed to have all four modifiers.

 Thank you for your help!

Version 1

See below:

modelH2a <- '

##DIRECT EFFECT

Y ~ c(M1YDirectA, M1YDirectB, M1YDirectC, M1YDirectD)*M1 + c(M2YDirectA, M2YDirectB, M2YDirectC, M2YDirectD)*M2 + c(M3YDirectA, M3YDirectB, M3YDirectC, M3YDirectD)*M3 + c(X1YDirectA, X1YDirectB, X1YDirectC, X1YDirectD)*X1 + c(X2YDirectA, X2YDirectB, X2YDirectC, X2YDirectD)*X2 + c(X3YDirectA, X3YDirectB, X3YDirectC, X3YDirectD)*X3 + c(X4YDirectA, X4YDirectB, X4YDirectC, X4YDirectD)*X4

##MEDIATORS

M1 ~ c(X1M1a, X1M1b, X1M1c, X1M1d)*X1 + c(X2M1a, X2M1b, X2M1c, X2M1d)*X2 + c(X3M1a, X3M1b, X3M1c, X3M1d)*X3 + c(X4M1a, X4M1b, X4M1c, X4M1d)*X4

M2 ~ c(X1M2a, X1M2b, X1M2c, X1M2d)*X1 + c(X2M2a, X2M2b, X2M2c, X2M2d)*X2 + c(X3M2a, X3M2b, X3M2c, X3M2d)*X3 + c(X4M2a, X4M2b, X4M2c, X4M2d)*X4

M3 ~ c(X1M3a, X1M3b, X1M3c, X1M3d)*X1 + c(X2M3a, X2M3b, X2M3c, X2M3d)*X2 + c(X3M3a, X3M3b, X3M3c, X3M3d)*X3 + c(X4M3a, X4M3b, X4M3c, X4M3d)*X4

##INDIRECT EFFECTS

###X1 -> Y via M1, M2, M3

abX1M1a := (X1M1a*M1YDirectA)

abX1M1b := (X1M1b*M1YDirectB)

abX1M1c := (X1M1c*M1YDirectC)

abX1M1d := (X1M1d*M1YDirectD)

abX1M2a := (X1M2a*M2YDirectA)

abX1M2b := (X1M2b*M2YDirectB)

abX1M2c := (X1M2c*M2YDirectC)

abX1M2d := (X1M2d*M2YDirectD)

abX1M3a := (X1M3a*M3YDirectA)

abX1M3b := (X1M3b*M3YDirectB)

abX1M3c := (X1M3c*M3YDirectC)

abX1M3d := (X1M3d*M3YDirectD)

##TOTAL EFFECTS

totalX1A := X1YDirectA + (X1M1a*M1YDirectA) + (X1M2a*M2YDirectA) + (X1M3a*M3YDirectA)

totalX1B := X1YDirectB + (X1M1b*M1YDirectB) + (X1M2b*M2YDirectB) + (X1M3b*M3YDirectB)

totalX1C := X1YDirectC + (X1M1c*M1YDirectC) + (X1M2c*M2YDirectC) + (X1M3c*M3YDirectC)

totalX1D := X1YDirectD + (X1M1d*M1YDirectD) + (X1M2d*M2YDirectD) + (X1M3d*M3YDirectD)'

Version 2

modelH2c <- '

##DIRECT EFFECT

Y ~ c(M1YDirectA, M1YDirectB, M1YDirectC, M1YDirectD)*M1 + c(M2YDirectA, M2YDirectB, M2YDirectC, M2YDirectD)*M2 + c(M3YDirectA, M3YDirectB, M3YDirectC, M3YDirectD)*M3 + c(X1YDirectA, X1YDirectB, X1YDirectC, X1YDirectD)*X1 + c(X2YDirectA, X2YDirectB, X2YDirectC, X2YDirectD)*X2 + c(X3YDirectA, X3YDirectB, X3YDirectC, X3YDirectD)*X3 + c(X4YDirectA, X4YDirectB, X4YDirectC, X4YDirectD)*X4

##MEDIATORS

M1 ~ c(X1M1a, X1M1b, X1M1c, X1M1d)*X1 + c(X2M1a, X2M1b, X2M1c, X2M1d)*X2 + c(X3M1a, X3M1b, X3M1c, X3M1d)*X3 + c(X4M1a, X4M1b, X4M1c, X4M1d)*X4

M2 ~ c(X1M2a, X1M2b, X1M2c, X1M2d)*X1 + c(X2M2a, X2M2b, X2M2c, X2M2d)*X2 + c(X3M2a, X3M2b, X3M2c, X3M2d)*X3 + c(X4M2a, X4M2b, X4M2c, X4M2d)*X4

M3 ~ c(X1M3a, X1M3b, X1M3c, X1M3d)*X1 + c(X2M3a, X2M3b, X2M3c, X2M3d)*X2 + c(X3M3a, X3M3b, X3M3c, X3M3d)*X3 + c(X4M3a, X4M3b, X4M3c, X4M3d)*X4

##INDIRECT EFFECTS

###X1 -> Y via M1, M2, M3

abM1 := (c(X1M1a, X1M1b, X1M1c, X1M1d)*X1) * (c(M1YDirectA, M1YDirectB, M1YDirectC, M1YDirectD)*M1)

abM2 := (c(X1M2a, X1M2b, X1M2c, X1M2d)*X1) * (c(M2YDirectA, M2YDirectB, M2YDirectC, M2YDirectD)*M2)

abM3 := (c(X1M3a, X1M3b, X1M3c, X1M3d)*X1) * (c(M3YDirectA, M3YDirectB, M3YDirectC, M3YDirectD)*M3)

##TOTAL EFFECTS

totalX1 := (c(X1YDirectA, X1YDirectB, X1YDirectC, X1YDirectD)*X1) + abM1 + abM2 + abM3'

[Mauricio Garnier-Villarreal]

my first guess is that you specified only one sample size. You need to specify a vector os sample sizes equal to the number of groups

For example

sample.nobs=c(100,100,100,100)

[Ash Gillis]

Thank you for your swift reply! Ok, this worked! However, I do have some additional questions

So, I tried to do this in simsem as well because it's my understanding I could obtain power through the simsem function 'getPower' once I run the simulation using the function 'sim'. 

However, I can't get the 'sim' function to run my simulation. I thought it was perhaps a similar issue to what you stated above. So I tried specifying a vector of sample sizes equal to the number of groups but no luck. 

Is it possible to obtain power in this way? Am I missing or misspecifying something in the sim argument? Or is there a way to use the result from the simulateData function in lavaan within the getPower function in simsem?

Thank you for any and all answers you can provide! 

Here is what I tried to do using simsem:

Sim1 <- sim(nRep=1000, model, n = c(250,250,250,250), lavaanfun = "sem")

Error in sim(nRep = 1000, fitH2a, n = c(250, 250, 250, 250), lavaanfun = "sem") : 

  When the number of replications argument 'nRep' is used, sample size 'n' must be a single value, a vector with length equal to 'nRep', or a vector with length equal to a multiple of 'nRep'

[Terrence Jorgensen]

The message indicates what the help page will tell you, which is that a vector for n= indicates the sample size to use in each replication.  With multigroup models, you need to specify a list (and if N varies across replications, specify a list of vectors).  This should work.

n = list(250,250,250,250)

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Mauricio Garnier-Villarreal]

simsem doent work for multiple gorups the same way with the lavaan syntax. You would need to use the matrix syntax to work that way.

With the lavaan syntaxyou need to cretae separta population data sets, like shown here https://github.com/simsem/simsem/blob/master/SupportingDocs/Examples/Version05/ex6/ex6.R

[Ash Gillis]

Ah ok. I tried Terrence's suggestion of using n = list(n,n,n,n). And I get the following error 

Error in names(Data) <- var.names : 

  'names' attribute [9] must be the same length as the vector [5]

Perhaps this is because of what Mauricio is saying about the current bug in simsem for multiproup models? Thank you Mauricio for sending this important info!

Just to make sure I have this straight Mauricio:

According to the example you sent, I would create separate population models for each group - do the parameters for this need to be prespecified or can they be freely estimated? They look prespecified in the example.

Then I also have the full multigroup model as I constructed in my first post and, in the sim function, use that original model in model =, and the dataframe with the results from the 4 population models as the 'rawdata'. Does that sound right?

Thank you so much!

Ash

[Mauricio Garnier-Villarreal]

If you only present the error message, without the syntax that produce it, its really hard for us to figure out what is the problem. We are guessing here what could be the problem, for example did you use n as the argument for sample size? because the argument for sample size sample.nobs.

Please share the full code that produce the errors

Also, its not a "bug" that simsem doesnt work the same for multiple groups with the lavaan syntax. Its a difference in features given how the user wants to use it

The main difference is that you have to create your full population first with simulatedata, like 10000 subjects per group first, then you pass this full population to simsem instead of a syntax. Then you have an analysis syntax, and a population data set for the simsem

[Ash Gillis]

Ah. My sincere apologies, I have been replying too quickly and not including essential info. I did use the n argument as opposed to sample.nobs. The sim() function sample size argument says n = like Terrance referred to above. When I use sample.nobs instead it provides an error saying I need to use n for sample size. 

Apologies for the term 'bug', the resource you sent refers to a bug but perhaps it relates to something different. I now understand how to follow the suggestion in the resource. 

Below is what I did. 

Simulation <- sim(nRep=800, n = list(200,200,200,200), fitH2a, data= PEB1, lavaanfun = "sem")

Error in names(Data) <- var.names : 

  'names' attribute [9] must be the same length as the vector [5]

Here is the model syntax as well, in case something in there is causing the issue. It is supposed to be a multigroup mediation model with three parallel mediators and three control variables (X2,X3,X4). I am eternally grateful for your help. Thank you.

fitH2a <- sem(modelH2a, data = PEB1, ordered= c("Y"), group = "Z")

[Ash Gillis]

  Looks like the syntax I just provided has been trimmed through gmail, so you may have to expand it to see. Thanks again

[Mauricio Garnier-Villarreal]

Also, I am not seing simulation parameters from your syntax. You need to provide population parameters for the data simulation

For simsem you need to simulate the data separately like

This would generate a large sample for 4 gorups

popData <- simulateData(popModel, sample.nobs = c(100000, 100000 , 100000 , 100000 ) )

And then run sim with the large population data, this way it would pick 200 random subjects from each

Output.a <- sim(1000, analyzeModel, n=list(200, 200), rawData = popData , group = "group", std.lv = TRUE, lavaanfun = "lavaan")

[Ash Gillis]

Yes yes, I was just providing syntax for the "analyzeModel" part

Thank you. I will return back later after constructing and running the popModel simulations and running sim() 

Thanks again

[Ash Gillis]

Hi there,

Thanks again for your previous help. I generated the population data for the four groups and ran sim with that data. This seemed to work.

I then tried to use this information to find power using varying sample sizes through getPower. My goal here is the find power for each parameter estimate for the actual sample sizes we have for each of the four groups. And to determine the sample sizes necessary for each parameter estimate at 80% power. 

However, I'm running into another issue. When I  use the getPower function, I'm provided with an object containing only one variable. Then when I pass that through findPower, I get the following error: Error in powerTable[, ivCol] : incorrect number of dimensions. I thought perhaps this was because the object I obtained through getPower was not a dataframe. So, I changed it to a dataframe and then received the following error when trying to use findPower:  Error in findPower(pow, "N", 0.8) : 

  Cannot find the specified target column. I have tried naming "N" as the variable name I get from the getPower dataframe (which is "V1") but I get the same error message. I appreciate any guidance you can provide. Thank you so much!

Here is my syntax for the population models, simulation, and power estimates

##Simulate population models for each group and combine into dataframe

popDatamodela <- simulateData(popModela, sample.nobs = 10000)

popDatamodelb <- simulateData(popModelb, sample.nobs = 10000)

popDatamodelc <- simulateData(popModelc, sample.nobs = 10000)

popDatamodeld <- simulateData(popModeld, sample.nobs = 10000)

popData <- rbind(popDatamodela, popDatamodelb, popDatamodelc, popDatamodeld)

popData <- data.frame(popData, group = rep(c(1, 2, 3, 4), each = 100000))

##estimate sample model from population data

Output.a <- sim(5, modelH2, n=list(200,200,200,200), rawData = popData, group = "group", std.lv = TRUE, lavaanfun = "sem") ##using nRep=5 in order to identify error, actual number will be much higher

##Obtain power

pow <- getPower(Output.a, alpha = 0.05, nVal=c(200,200,200,200), pmMCAR = c(0, 0.1, 0.2)) ###The 'pow' object/dataframe contains only one variable

findPower(pow, "N", 0.80) 

>>Error in powerTable[, ivCol] : incorrect number of dimensions

pow <- data.frame(pow)

findPower(pow, "N", 0.80)

>>Error in findPower(pow, "N", 0.8) : Cannot find the specified target column

[Mauricio Garnier-Villarreal]

To use the power analysis functionality to have to estimate the models with varying sample sizes. So the functions can calculate the power at different sample sizes

Output.a <- sim(NULL, modelH2, n=list(50:500, 50:500 , 50:500 , 50:500 ), rawData = popData, group = "group", std.lv = TRUE, lavaanfun = "sem")

[Ash Gillis]

Thank you so much! This worked for me.

I just have one last question.

I see from my result of getPower that some power estimates are equal to 1 for some parameters. Is there any reason for why this may be the case? 

I am calculating power from small samples as you'll see from nVal= so it seems to me unlikely that power is so high it would be equal to .999

Output.2 <- sim(51, modelH2, n=list(50:100,50:100,50:100,50:100), rawData = popData, group = 'group', std.lv = TRUE, lavaanfun = "sem")

summary(Output.2)

##Find power for mediation models at actual sample sizes for each group

pow <- getPower(Output.2, alpha = 0.05, nVal=c(76, 50, 44, 59))

[Mauricio Garnier-Villarreal]

this should be because these parameters have strong effect sizes, like high regression slopes

[Alex Schoemann]

A couple of thoughts about your approach here.

I would make sure the simulation includes all sample sizes you're planning to compute power for. It looks like one of your group sizes is 44, but the smallest sample size in a group you're looking at is 50. Ideally the smallest sample size used in the simulation would be below the smallest N of interest.

In all your simulations you're assuming equal group sizes for each replication (N in all groups is either 50, 51, etc.)*

You're only using 51 replications. For stable estimates of power you need many more replications, ideally 1000+ or even 5000+. You could change the n option to be: list(rep(50:100,100),rep(50:100,100),rep(50:100,100),rep(50:100,100)) this will give you 5100 replications.

* The way getPower is currently built it will only use the overall N from the model to compute power in a situation where N varies across replications. This isn't ideal for a multiple group model (and it's something we need to fix) since you can arrive at the same total N with different N per group and that might have an impact on power. I think for now your approach is a good one, but be aware that your power analysis has an assumption that doesn't match your data.

Alex

[Ash Gillis]

Hi Alex,

Thanks for your comments. 

I will go ahead and make the smallest sample size in the simulation below the smallest N of interest, thanks for catching that.

I'll make a note about the power analysis' assumption of equal group sizes for each replication.

Also, yes, I used the 51 replications just so I could identify errors I was receiving previously without running thousands of iterations. I plan to use at least 5000. 

Thanks again!

And special thanks to Mauricio who helped me throughout this process! I am eternally grateful

Ash

[Gb]

Hi!

I do have a similar issue in running my simulation to detect the sample size for a model I want to test, but I got this error message: Error in sim(NULL, model = analysis.mod1, n = list(500, 500), rawData = popData, : Check the list object specified in the 'rawData' argument; list must either contain matrices or data frames In addition: Warning message: In sim(NULL, model = analysis.mod1, n = list(500, 500), rawData = popData, : The sample size argument 'n' is suppressed when 'rawData' is specified.

I am sure about how to write the analysis mod. Plus, not sure if I have to run the simulation for the 4 types of invariances (noninvariance, configural, weak and strong). Can you help me out?

#simulate covariates
gender <- rep(c(0,1), 800)
age <- rep (18:80, 25)
fs <-rep(c(0,1), 800)
ed <- rep(1:7, 200)
nationality <- rep(c(0,1), 800)
income <- rep(1:10, 160)


#run pop model with 6 covariates included  
#modpopa

popModela <- "
#measurament model
sti =~ 0.70*s1 + 0.70*s2 + 0.70*s3
pc =~ 0.70*p1 + 0.70*p2 + 0.70*p3
ic =~ 0.70*i1 + 0.70*i2 + 0.70*i3
att =~ 0.70*a1 + 0.70*a2 + 0.70*a3
sn =~ 0.70*sn1 + 0.70*sn2 + 0.70*sn3
pbc =~ 0.70*pb1 + 0.70*pb2 + 0.70*pb3 + 0.70*pb4
int =~ 0.70*in1 + 0.70*in2 + 0.70*in3
pb =~ 0.70*pb1


pb~~ 1*pb                
sti ~~ 1*sti
pc ~~ 1*pc
ic ~~ 1*ic
att  ~~ 1*att
sn ~~ 1*sn
pbc ~~ 1*pbc
int ~~ 1*int
 
sti  ~~ 0.5*pc
sti ~~ 0.5*ic
sti ~~ 0.5*att
sti ~~ 0.5*sn
sti ~~ 0.5*pbc
sti ~~ 0.5*int
sti ~~ 0.5*pb
pc ~~ 0.5*ic
pc ~~ 0.5*att
pc ~~ 0.5*sn
pc ~~ 0.5*pbc
pc ~~ 0.5*int
ic ~~ 0.5*att
ic ~~ 0.5*sn
ic ~~ 0.5*pbc
ic ~~ 0.5*int
ic ~~ 0.5*pb
att ~~ 0.5*sn
att ~~ 0.5*pbc
att ~~ 0.5*int
att ~~ 0.5*pb
sn ~~ 0.5*pbc
sn ~~ 0.5*int
sn ~~ 0.5*pb
pbc ~~ 0.5*int
pbc ~~ 0.5*pb
int ~~ 0.5*pb

s1 ~~ 0.49*s1
s2 ~~ 0.49*s2
s3 ~~ 0.49*s3
p1 ~~ 0.49*p1
p2 ~~ 0.49*p2
p3 ~~ 0.49*p3
i1 ~~ 0.49*i1
i2 ~~ 0.49*i2
i3 ~~ 0.49*i3
a1 ~~ 0.49*a1
a2 ~~ 0.49*a2
a3 ~~ 0.49*a3
sn1 ~~ 0.49*sn1
sn2 ~~ 0.49*sn2
sn3 ~~ 0.49*sn3
pbc1 ~~ 0.49*pbc1
pbc2 ~~ 0.49*pbc2
pbc3 ~~ 0.49*pbc3
pbc4  ~~ 0.49*pbc4
in1~~ 0.49*in1
in2~~ 0.49*in2
in3~~ 0.49*in3
pb1 ~~ 1*pb1

#structural model


att ~ -0.50*sti  + 0.50*pc + 0.50*ic
int ~ 0.50*att + 0.50*sn + 0.50*pbc + 0.50*pb
"
#popmodb
popModelb <- "
#measurament model
sti =~ 0.70*s1 + 0.70*s2 + 0.70*s3
pc =~ 0.70*p1 + 0.70*p2 + 0.70*p3
ic =~ 0.70*i1 + 0.70*i2 + 0.70*i3
att =~ 0.70*a1 + 0.70*a2 + 0.70*a3
sn =~ 0.70*sn1 + 0.70*sn2 + 0.70*sn3
pbc =~ 0.70*pb1 + 0.70*pb2 + 0.70*pb3 + 0.70*pb4
int =~ 0.70*in1 + 0.70*in2 + 0.70*in3
pb =~ 0.70*pb1


pb~~ 1*pb                
sti ~~ 1*sti
pc ~~ 1*pc
ic ~~ 1*ic
att  ~~ 1*att
sn ~~ 1*sn
pbc ~~ 1*pbc
int ~~ 1*int
 
sti  ~~ 0.5*pc
sti ~~ 0.5*ic
sti ~~ 0.5*att
sti ~~ 0.5*sn
sti ~~ 0.5*pbc
sti ~~ 0.5*int
sti ~~ 0.5*pb
pc ~~ 0.5*ic
pc ~~ 0.5*att
pc ~~ 0.5*sn
pc ~~ 0.5*pbc
pc ~~ 0.5*int
ic ~~ 0.5*att
ic ~~ 0.5*sn
ic ~~ 0.5*pbc
ic ~~ 0.5*int
ic ~~ 0.5*pb
att ~~ 0.5*sn
att ~~ 0.5*pbc
att ~~ 0.5*int
att ~~ 0.5*pb
sn ~~ 0.5*pbc
sn ~~ 0.5*int
sn ~~ 0.5*pb
pbc ~~ 0.5*int
pbc ~~ 0.5*pb
int ~~ 0.5*pb

s1 ~~ 0.49*s1
s2 ~~ 0.49*s2
s3 ~~ 0.49*s3
p1 ~~ 0.49*p1
p2 ~~ 0.49*p2
p3 ~~ 0.49*p3
i1 ~~ 0.49*i1
i2 ~~ 0.49*i2
i3 ~~ 0.49*i3
a1 ~~ 0.49*a1
a2 ~~ 0.49*a2
a3 ~~ 0.49*a3
sn1 ~~ 0.49*sn1
sn2 ~~ 0.49*sn2
sn3 ~~ 0.49*sn3
pbc1 ~~ 0.49*pbc1
pbc2 ~~ 0.49*pbc2
pbc3 ~~ 0.49*pbc3
pbc4  ~~ 0.49*pbc4
in1~~ 0.49*in1
in2~~ 0.49*in2
in3~~ 0.49*in3
pb1 ~~ 1*pb1

#structural model


att ~ -0.50*sti  + 0.50*pc + 0.50*ic
int ~ 0.50*att + 0.50*sn + 0.50*pbc + 0.50*pb
"

#simulate population model for each group and combine into dataframe

popDatamodela <- simulateData(popModela, sample.nobs = 10000)
popDatamodelb <- simulateData(popModelb, sample.nobs = 10000)

popData <- rbind(popDatamodela, popDatamodelb)
popData <- data.frame(popData, group = rep(c(1, 2), each = 100000))




analysis.mod1  <- "
# define latent variables
sti =~ s1 + s2 + s3
pc =~ p1 + p2 + p3
ic =~ i1 + i2 + i3
att =~ a1 + a2 + a3
sn =~ sn1 + sn2 + sn3
pbc =~ pb1 + pb2 + pb3 + pb4
int =~ in1 + in2 + in3
pb =~ pb1


# specify structural model

att ~ c(b8, z8)*sti+ c(b10,z10)*pc + c(b11, z11)*ic
int ~ c(b2, z2)*att + c(b4, z4)*sn + c(b6, z6)*pbc + c(b7, z7)*pb + c(b1, z1)*sti + c(b12, z12)*pc + c(b3, z3)*ic

# indirect effects group a
st.att.int := b8*b2
pc.att.int := b10*b2
ic.att.int := b11*b2

# indirect effects group b
st.att.int := z8*z2
pc.att.int := z10*z2
ic.att.int := z11*z2

# total effects group a
                 
tt.st := b1 + b8*b2
tt.pc := b12 +b10*b2
tt.ic := b3 + b11*b2

# total effects group b
                 
tt.st := z1 + z8*z2
tt.pc := z12 +z10*z2
tt.ic := z3 + z11*z2
"
med.n1 <- sim(NULL,  model = analysis.mod1, n=list(500,500), rawData = popData, group = "group", std.lv = TRUE, estimator = "MLR", lavaanfun = "sem")
summary(med.n1)
powTable2 <- getPower(med.n1)
findPower(powTable2, "N", 0.85)

[Gb]

Sorry I forgot to add the covariates on intention that it would be like that:

0.5*gender + 0.5*age + 0.5*fs +0.5*nationality + 0.5*income (for each of the pop model)

c(b12, z12)*gender + c(b13, z13)*age + c(b14,z14)*fs +(b15, z15)*nationality + c(b16,z16)*income (this for the analysis mod). 

Not sure if I imputed the right way the covariates (probably not)

[Terrence Jorgensen]

Not sure if I imputed the right way the covariates (probably not)

These vignettes show how to incorporate exogenous covariates:

https://github.com/simsem/simsem/blob/master/SupportingDocs/Examples/Version05/exMultipleFormat/covData.R

https://github.com/simsem/simsem/blob/master/SupportingDocs/Examples/Version05/exMultipleFormat/moderatedMediation.R

[HG Wells]

Hi everyone, I am trying to run a similar example as per the original post, but with two groups; I encouter some error messages.

This works:

# Simulation
sim_out <- sim(5000,
               test_model,
               n = list(100, 100),

               rawData = popData,
               group = "group",
               std.lv = TRUE,

               lavaanfun = "sem",
               seed = 1234)

But when I try to implement Alex's suggestion: 

You're only using 51 replications. For stable estimates of power you need many more replications, ideally 1000+ or even 5000+. You could change the n option to be: list(rep(50:100,100),rep(50:100,100),rep(50:100,100),rep(50:100,100)) this will give you 5100 replications.

I run this and get the error message, below:

sim_out <- sim(NULL,
               test_model,
               n = list(rep(50:100, 100), rep(50:100, 100)),

               rawData = popData,
               group = "group",
               std.lv = TRUE,

               lavaanfun = "sem",
               seed = 1234)

I get the following error: 

Error in sim(NULL, test_model, n = list(rep(50:100, 100), rep(50:100, : Check the list object specified in the 'rawData' argument; list must either contain matrices or data frames In addition: Warning message: In sim(NULL, test_model, n = list(rep(50:100, 100), rep(50:100, : The sample size argument 'n' is suppressed when 'rawData' is specified

What do I do wrong? How I should I tweak the syntax to obtain those 5,100 replications?

Thank you in advance for any help you may be able to offer. 

Emanuele