Using lavaan within a function

[Jarrod]

Hi, I'm attempting to create a function that would allow me to run several iterations of the same model, and then table the parameter estimates of those models; however, I am running into an issue that only appears within the function (not when I use the code outside of the function). To start with a basic example, I'm just attempting to run a univariate regression model within this function:

lavaan.uv.reg <- function(outcome, predictor, RegDat) {
  reg.model <- '
                outcome ~ predictor
                outcome ~~ outcome
                outcome ~ 1
                '
  reg.run <- sem(model = reg.model, data = RegDat)
  est <- reg.run@Fit@est[1]
  se <- reg.run@Fit@se[1]
  z <- est/se
  p <- 2*pnorm(-abs(reg.run.z))
  reg.run.output <- cbind(est, se, z, p)
  return(reg.run.output)
}
TableA1 <- lavaan.uv.reg(outcomeA, predictor1, mydata)
TableA2 <- lavaan.uv.reg(outcomeA, predictor2, mydata)
TabledResults <- rbind(TableA1, TableA2)


When I run this code, I get the following error, indicating that my variables are not in my dataset:

TableA1 <- lavaan.uv.reg(outcomeA, predictor1, mydata)
 Error in lav_data_full(data = data, group = group, group.label = group.label,  :
  lavaan ERROR: missing observed variables in dataset: outcome predictor

I know this code indicates that 1+ variables are not in the dataset, but the code runs fine if I extract it from the function and run it on its own. I'm assuming this issue has something to do with trying to using this code within a function (maybe the way functions handle code within ' ' ?), but I haven't been able to specify the issue or whether its possible to incorporate this code into a new, user-defined, function.

Any help would be greatly appreciated.
Thank you,
Jarrod

sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.6 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base    

other attached packages:
 [1] knitr_1.14            plyr_1.8.4            GPArotation_2014.11-1 Rcsdp_0.1.55          psych_1.6.9           OpenMx_2.6.9        
 [7] Rcpp_0.12.7           Matrix_1.2-7.1        MASS_7.3-45           digest_0.6.10         lavaan_0.5-22       

loaded via a namespace (and not attached):
[1] quadprog_1.5-5  lattice_0.20-34 grid_3.3.1      stats4_3.3.1    pbivnorm_0.6.0  tools_3.3.1     foreign_0.8-67  rsconnect_0.4.3 mnormt_1.5-4  

[Jeremy Miles]

I think the variable names inside quotes are ignored, and don't change.

Instead of:

  reg.model <- '
                outcome ~ predictor
                outcome ~~ outcome
                outcome ~ 1
                '

 

Use something like:

  reg.model <-
                paste0(outcome,  " ~ ", predictor ' "\n",
                outcome, " ~~ ", outcome, "\n"
                outcome, " ~ 1")

(Note: untested code.)

Hope that helps. If it doesn't, a reproducible example would be easier to fix.

Jeremy               

 

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at https://groups.google.com/group/lavaan.
For more options, visit https://groups.google.com/d/optout.

[Terrence Jorgensen]

I think the variable names inside quotes are ignored, and don't change.

Jeremy is right, but his solution would only work if the first two arguments of your functions were character strings, not vectors of values.  It is not really clear which you need from your code.  In fact, the first two arguments seem unnecessary.  If the variables in the data.frame you pass to the RegDat argument (i.e., mydata) are always going to be named "outcome" and "predictor", then the only argument your function needs is RegDat.  The main point is that lavaan will read your model syntax as is, look at the variable names it contains, and look for columns in its data argument that correspond to the (non-latent) variable names.  

If the only aspects that change between iterations are the values of the variables, then your function could instead include arguments about data generation, and you could create RegDat (a data.frame with columns "outcome" and "predictor") within the function before fitting the (constant) model to the data.  If you search the topics in this group for "simulate", you might find some helpful simulation examples in previous posts.

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

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

UvA web page: http://www.uva.nl/profile/t.d.jorgensen

[kma...@aol.com]

Jarod,
  I will offer two further thoughts to complement the strategies offered by Jeremy and Terrence.  (Careful, I think that there may be a single quote in place of a comma in Jeremy's example).

(1)  When using the paste() strategy suggested by Jeremy, I like to enclose expressions to be evaluated in the eval() function.  This is not necessary and some may consider it excessive but I include it as a form of defensive programming to avoid use/mention errors in the code.

(2) Remember that it is possible to feed a parameter table instead of a model to one of lavaan's model fitting functions.

AA <- '[model syntax]'
BB <- lavaanify(AA)
CC <- lavaan(BB, ...)

Sometimes it is more convenient, and produces more readable and extensible code, to edit a parameter table instead of model syntax.  You can adopt a hybrid strategy of generating a parameter table with lavaanify() applied to template model syntax, and then modifying the parameter table (including variable names).  Alternatively, you can directly construct a parameter table.

DD <- data.frame(
  id=NA,
  lhs=NA,
  op=NA,
  rhs=NA,
  user=NA,
  group=NA,
  free=NA,
  ustart=NA,
  exo=NA,
  label=NA,
  plabel=NA)
DD
newLine <- c(1,'y','~','x',1,1,1,NA,0,'','.p1')
DD[1,] <- newLine
DD

Your choice of strategy might depend in part on both the complexity of the models and the extent to which you anticipate extending and reusing the resulting function.  If you can define n-1 of your models as predictable modifications of the nth model, you might consider making the model syntax for the nth model a parameter to your function, lavaanifying it within the function, and then generating the other models by modifying the resulting parameter table.

Keith Markus

[Yves Rosseel]

> Alternatively, you can
> directly construct a parameter table.
>
> DD <- data.frame(
> id=NA,
> lhs=NA,
> op=NA,
> rhs=NA,
> user=NA,
> group=NA,
> free=NA,
> ustart=NA,
> exo=NA,
> label=NA,
> plabel=NA)

In recent versions of lavaan (> 0.5-21), only a few columns are required
if you provide a parameter table. Essential columns are "lhs", "op" and
"rhs". In this case, we assume everything is free. If you add a 'free'
column, non-zero numbers (often just 1) indicate free values, everything
else is fixed to zero. Except if you add a 'ustart' column, which is
useful to provide non-zero but fixed parameters.

As an example:

example(cfa)
PT <- parTable(fit)[,c("lhs","op","rhs","free","ustart")]
PT$free[ PT$free > 0] <- 1
PT

fit2 <- sem(PT, data = HolzingerSwineford1939)
fit2

lavaan (0.5-23.1020) converged normally after 35 iterations

Number of observations 301

Estimator ML
Minimum Function Test Statistic 85.306
Degrees of freedom 24
P-value (Chi-square) 0.000


Yves.

[kma...@aol.com]

Yves' welcome simplification makes this clarification largely moot.  However, if it can help save others time debugging, I probably should have mentioned that ustart needs to be encoded as numeric (NA values may default to character) and group and exo need to be encoded as integer (integer values may default to numeric).

DD$ustart <- as.numeric(DD$ustart)
DD$exo <- as.integer(DD$exo)
DD$group <- as.integer(DD$group)

Keith Markus

[Yves Rosseel]

Good point. I took care of this automatically in dev 0.5-23.1021

Yves.

[Jarrod]

Hi All, Thanks for your help on this. Jeremy's fix worked, and in regard to Terrence's comment, I used many different IVs and DVs. Was just trying to develop a quick bit of code (starting simple) to extract parameters that I could table. This exercise was, in part, to help me get more familiar with R and also to help me learn more about lavaan. 

The final, functional code that I plan to adapt is pasted below, though I'll also take a closer look at the functions mentioned by Keith (thanks!).

lavaan.uv.reg <- function(outcome, predictor, RegData) {
  reg.model <-  paste0(outcome,  " ~ ", predictor, "\n",
                outcome, " ~~ ", outcome, "\n",
                outcome, " ~ 1")
  reg.run <- sem(model = reg.model, data = RegData)

  est <- reg.run@Fit@est[1]
  se <- reg.run@Fit@se[1]
  z <- est/se

  p <- 2*pnorm(-abs(z))