Parameter constraints involving power relations

[serena defina]

Good afternoon, 

I am constructing a CLPM and would like to specify equality constraints that account for different time intervals between occasions, i.e., i would like the stability over any number of intervals to be a function of the stability across one interval raised to the power of the number of intervals that have elapsed

So, for example if 1 year passed between t1 and t2 but one and a half years passed between t2 and t3, I would like these relations to be expressed as:

b2 = b1^1.5 

where b1 is the relation between t1 and t2 and b2 is the relation between t2 and t3.

I can specify this successfully in lavaan (0.6-16) using integer exponents, but when a float is imputed (no matter the precision) I get:

`In lavaan::lavaan(model = formula(dep_temp, cmr_temp, mod), data = data,  :
  lavaan WARNING:
    Model estimation FAILED! Returning starting values.`

I have tried using `**` instead of `^` (weak attempt, i know). Any idea how I could solve or work around this?

Here is a copy of an the formula I use (note the constraints section at the end) 

# Unit effects
eta_dep =~ 1*dep1 + 1*dep2 + 1*dep3 + 1*dep4 + 1*dep5 + 1*dep6
eta_cmr =~ 1*cmr1 + 1*cmr2 + 1*cmr3 + 1*cmr4 + 1*cmr5 + 1*cmr6
# Impulses
u_dep1 =~ dep1
dep1 ~~ 0*dep1
u_dep2 =~ dep2
dep2 ~~ 0*dep2
u_dep3 =~ dep3
dep3 ~~ 0*dep3
u_dep4 =~ dep4
dep4 ~~ 0*dep4
u_dep5 =~ dep5
dep5 ~~ 0*dep5
u_dep6 =~ dep6
dep6 ~~ 0*dep6
u_cmr1 =~ cmr1
cmr1 ~~ 0*cmr1
u_cmr2 =~ cmr2
cmr2 ~~ 0*cmr2
u_cmr3 =~ cmr3
cmr3 ~~ 0*cmr3
u_cmr4 =~ cmr4
cmr4 ~~ 0*cmr4
u_cmr5 =~ cmr5
cmr5 ~~ 0*cmr5
u_cmr6 =~ cmr6
cmr6 ~~ 0*cmr6
# Regressions
dep6 ~ AR_dep5*dep5 + CL_dep5*cmr5 + maAR_dep5*u_dep5 + maCL_dep5*u_cmr5
dep5 ~ AR_dep4*dep4 + CL_dep4*cmr4 + maAR_dep4*u_dep4 + maCL_dep4*u_cmr4
dep4 ~ AR_dep3*dep3 + CL_dep3*cmr3 + maAR_dep3*u_dep3 + maCL_dep3*u_cmr3
dep3 ~ AR_dep2*dep2 + CL_dep2*cmr2 + maAR_dep2*u_dep2 + maCL_dep2*u_cmr2
dep2 ~ AR_dep1*dep1 + CL_dep1*cmr1 + maAR_dep1*u_dep1 + maCL_dep1*u_cmr1
cmr6 ~ CL_cmr5*dep5 + AR_cmr5*cmr5 + maCL_cmr5*u_dep5 + maAR_cmr5*u_cmr5
cmr5 ~ CL_cmr4*dep4 + AR_cmr4*cmr4 + maCL_cmr4*u_dep4 + maAR_cmr4*u_cmr4
cmr4 ~ CL_cmr3*dep3 + AR_cmr3*cmr3 + maCL_cmr3*u_dep3 + maAR_cmr3*u_cmr3
cmr3 ~ CL_cmr2*dep2 + AR_cmr2*cmr2 + maCL_cmr2*u_dep2 + maAR_cmr2*u_cmr2
cmr2 ~ CL_cmr1*dep1 + AR_cmr1*cmr1 + maCL_cmr1*u_dep1 + maAR_cmr1*u_cmr1
# Comovement
u_dep1 ~~ comv1*u_cmr1
u_dep2 ~~ comv2*u_cmr2
u_dep3 ~~ comv3*u_cmr3
u_dep4 ~~ comv4*u_cmr4
u_dep5 ~~ comv5*u_cmr5
u_dep6 ~~ comv6*u_cmr6
# Restrictions
u_dep1 ~~ 0*eta_dep + 0*eta_cmr + 0*u_dep2 + 0*u_dep3 + 0*u_dep4 + 0*u_dep5 + 0*u_dep6 + 0*u_cmr2 + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_dep2 ~~ 0*eta_dep + 0*eta_cmr + 0*u_dep3 + 0*u_dep4 + 0*u_dep5 + 0*u_dep6 + 0*u_cmr1 + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_dep3 ~~ 0*eta_dep + 0*eta_cmr + 0*u_dep4 + 0*u_dep5 + 0*u_dep6 + 0*u_cmr1 + 0*u_cmr2 + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_dep4 ~~ 0*eta_dep + 0*eta_cmr + 0*u_dep5 + 0*u_dep6 + 0*u_cmr1 + 0*u_cmr2 + 0*u_cmr3 + 0*u_cmr5 + 0*u_cmr6
u_dep5 ~~ 0*eta_dep + 0*eta_cmr + 0*u_dep6 + 0*u_cmr1 + 0*u_cmr2 + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr6
u_dep6 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr1 + 0*u_cmr2 + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr5
u_cmr1 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr2 + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_cmr2 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr3 + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_cmr3 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr4 + 0*u_cmr5 + 0*u_cmr6
u_cmr4 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr5 + 0*u_cmr6
u_cmr5 ~~ 0*eta_dep + 0*eta_cmr + 0*u_cmr6
u_cmr6 ~~ 0*eta_dep + 0*eta_cmr
# Constraints
AR_dep2 == AR_dep1^0.5
AR_cmr2 == AR_cmr1^1
CL_dep2 == CL_dep1^0.5
CL_cmr2 == CL_cmr1^1
AR_dep3 == AR_dep1^1.3
AR_cmr3 == AR_cmr1^0.8
CL_dep3 == CL_dep1^1.3
CL_cmr3 == CL_cmr1^0.8
AR_dep4 == AR_dep1^0.5
AR_cmr4 == AR_cmr1^1.2
CL_dep4 == CL_dep1^0.5
CL_cmr4 == CL_cmr1^1.2
AR_dep5 == AR_dep1^2.7
AR_cmr5 == AR_cmr1^3.3
CL_dep5 == CL_dep1^2.7
CL_cmr5 == CL_cmr1^3.3

[Terrence Jorgensen]

I can specify this successfully in lavaan (0.6-16) using integer exponents, but not when a float is imputed (no matter the precision)

Try installing the latest version:  remotes::install_github("yrosseel/lavaan")

I can use noninteger exponents in this simpler example:

> mod <- "x3 ~ b1*x1 + b2*x2 ; b2 == b1^1.5"
> fit <- sem(mod, data = HolzingerSwineford1939)
> summary(fit)
lavaan 0.6.17.1902 ended normally after 14 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         3

  Number of observations                           301

Model Test User Model:
                                                     
  Test statistic                                 0.002
  Degrees of freedom                                 1
  P-value (Chi-square)                           0.967

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  x3 ~                                                
    x1        (b1)    0.363    0.032   11.485    0.000
    x2        (b2)    0.218    0.029    7.657    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .x3                0.966    0.079   12.268    0.000

Constraints:
                                               |Slack|
    b2 - (b1^1.5)                                0.000

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

[serena defina]

Thanks Terrence! 

Unfortunately no better luck with the latest version (0.6-17.1966) either, even with a minimal (clp) model. I suppose then this has to do with the number of constraints that are imposed...?

best, 

Serena 

[Shu Fai Cheung]

May you post the call you used to fit the model (e.g., the call to sem())? It may help figuring out the causing of the problem.

-- Shu Fai

--
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 view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/583746c4-0e8f-493a-8580-f28638c79776n%40googlegroups.com.

[serena defina]

Sure, I am calling lavaan() actually ( but encountered the same error with sem() ) 

m <- lavaan::lavaan( riclpm_formula(meas_time = dc[3:4], stationarity = stationarity), # output: string specifying formula
                      data = d,
                      se = 'robust',
                      missing = 'fiml', # group = group,
                      meanstructure = TRUE,
                      int.ov.free = TRUE
                      )

[Shu Fai Cheung]

I have never used riclpm_formula(). What package or source does it come from? Is it the function used to generate the model syntax you posted?

-- Shu Fai

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/3d98be71-f964-4b1a-b3ab-ffcab8d77c03n%40googlegroups.com.

[serena defina]

Sorry yes that is user defined but it simply outputs the model syntax, based on the RI_CLPM structure described here https://jeroendmulder.github.io/RI-CLPM/lavaan.html for example 

I pasted the function here for reference (highlighted in red the portion where I set the equality constraints I want to weight 

So for example if you now call 

f = riclpm_formula('x', 'y', meas_time=list(c(3, 4.5, 5.2), c(2.9, 4.4, 6)))

it should print out an example where x was measured at 3, 4.5, and 5.2 years and y was measured at 2.9, 4.4, and 6 years. 

riclpm_formula <- function(var1='dep', var2='cmr',

                                                n_ocs=NULL, meas_time=list(c(), c()),
                                                stationarity=TRUE,

                                                verbose=TRUE) {

 
  # Check input ----------------------------------------------------------------
  if (is.null(n_ocs) & (length(meas_time[[1]])==0 | length(meas_time[[1]])==0)) {
    stop('Provide total number of occasions or time points of measurement!') }
 
  if (length(meas_time[[1]]) != length(meas_time[[2]])) {
    stop('You need the same number of time points for each of the two variables.') }
 
  if (!is.null(n_ocs) & !length(meas_time[[1]]) %in% c(0, n_ocs)) { # contraddictory info
    stop('Number of occations and measurement times provided do not agree.') }
 
  # ----------------------------------------------------------------------------
  # How many occasions
  if (length(meas_time[[1]]) > 0) {
    temp_var1 = meas_time[[1]]
    temp_var2 = meas_time[[2]]
  } else {
    temp_var1 = temp_var2 = 1:n_ocs
  }
 
  if (is.null(n_ocs)) { n_ocs = length(temp_var1) }
 
  # Create between components (random intercepts)
  random_intercepts = paste0('# Random intercepts\n',
                             'ri_',var1,' =~ ', paste0('1*',var1,1:n_ocs, collapse=' + '), '\n',
                             'ri_',var2,' =~ ', paste0('1*',var2,1:n_ocs, collapse=' + '), '\n')
 
  # Create within-components (or within-person centered variables)
  # Note: factor loadings are set to 1 in non-stationary model and freely estimated in stationary models
  if (stationarity) { l = 'NA*' } else { l = '1*' }
 
  impulses = paste0('# Impulses\n',
                    paste0('w_',var1,1:n_ocs,' =~ ',l,var1,1:n_ocs, collapse='\n'), '\n',
                    paste0('w_',var2,1:n_ocs,' =~ ',l,var2,1:n_ocs, collapse='\n'), '\n')
 
  # Estimate lagged effects between within-person centered variables
  regressions = '# Regressions\n'
  for (i in 2:n_ocs) {
    # if (n_strata > 1) { ... ar_1 <- paste0('c(', paste0(rep(paste0('AR_',var1,i-1),strata), collapse=', '), ')') #define  strata

    ar_1 <- paste0('AR_',var1,i-1)
    ar_2 <- paste0('AR_',var2,i-1)
    cl_1 <- paste0('CL_',var1,i-1)
    cl_2 <- paste0('CL_',var2,i-1)
   
    regressions = paste0(regressions,
      'w_',var1,i,' ~ ',ar_1,' * w_',var1,i-1,' + ',cl_1,' * w_',var2,i-1,'\n',
      'w_',var2,i,' ~ ',ar_2,' * w_',var2,i-1,' + ',cl_2,' * w_',var1,i-1,'\n')
  }
 
  # Estimate covariance between within-person centered variables
  covariances = paste0('# Covariances\n',
                       # Estimate correlation between within-components at first wave
                       'w_',var1,'1 ~~ cor1 * w_',var2,'1\n',

                       # Estimate residual covariances (between residuals of within-person vars, i.e., innovations)
                       paste0('w_',var1,2:n_ocs,' ~~ rcov',2:n_ocs,' * w_',var2,2:n_ocs, collapse='\n'), '\n',

                       # Estimate covariance of the random intercepts
                       'ri_',var1,' ~~ covRI * ri_',var2,'\n')
 
  if (stationarity) { w1 = '1*' } else { w1 = '' }
 
  variances = paste0('# Variances\n',
                     # Set variances of within-components at first wave to 1
                     paste0('w_',var1,'1 ~~ ',w1,'w_',var1,'1\n', 'w_',var2,'1 ~~ ',w1,'w_',var2,'1\n'),

                     # Estimate (and label) residual variances of within-person centered variables
                     paste0('w_',var1,2:n_ocs,' ~~ rvar_',var1,2:n_ocs,'*w_',var1,2:n_ocs, collapse='\n'),'\n',
                     paste0('w_',var2,2:n_ocs,' ~~ rvar_',var2,2:n_ocs,'*w_',var2,2:n_ocs, collapse='\n'),'\n',

                     # Estimate variance of random intercepts
                     'ri_',var1,' ~~ ri_',var1,'\nri_',var2,' ~~ ri_',var2, '\n')
 
  constraints = ar_var1_con = ar_var2_con = cl_var1_con = cl_var2_con = cor_con = rvar_con = ''
 
  if (stationarity) {
    constraints = '# Constraints\n'
     
    # paste0('# Constrain grand means over time\n',
    # paste0(var1,1:n_ocs, collapse=' + '), ' ~ mean_',var1,'*1\n',
    # paste0(var2,1:n_ocs, collapse=' + '), ' ~ mean_',var2,'*1\n')
   
    # Weight for unequal temporal gaps
    rel='^' # '*'
    i = 1:(n_ocs-2) # iterator

    prec=2 # precision
    # NOTE: round to nearest integer for the moment 
    ar_var1_con <- paste0('AR_',var1,i,rel,round(temp_var1[i+1]-temp_var1[i], prec), ' == AR_',var1,i+1,rel,round(temp_var1[i+2]-temp_var1[i+1], prec),

                          collapse='\n')
    ar_var2_con <- paste0('AR_',var2,i,rel,round(temp_var2[i+1]-temp_var2[i], prec), ' == AR_',var2,i+1,rel,round(temp_var2[i+2]-temp_var2[i+1], prec),

                          collapse='\n')
    cl_var1_con <- paste0('CL_',var1,i,rel,round(temp_var1[i+1]-temp_var2[i], prec), ' == CL_',var1,i+1,rel,round(temp_var1[i+2]-temp_var2[i+1], prec),

                          collapse='\n')
    cl_var2_con <- paste0('CL_',var2,i,rel,round(temp_var2[i+1]-temp_var1[i], prec), ' == CL_',var2,i+1,rel,round(temp_var2[i+2]-temp_var1[i+1], prec),

                          collapse='\n')
   
    # Compute correlations between the within-components themselves at each wave
    i = 1:(n_ocs-1) # new iterator
    cor_con <- paste0('cor',2:n_ocs, ' := AR_',var1,i,'*CL_',var2,i,' + ', 'CL_',var1,i,'*AR_',var2,i,' + ',  

                                                                  'AR_',var1,i,'*AR_',var2,i,'*cor',i,' + ', 'CL_',var1,i,'*CL_',var2,i,'*cor',i,' + rcov',2:n_ocs, collapse='\n')
   
    # Constrain residual variances of within-components such that the total variance of each
    # within-component equals 1
    rvar_con <- paste0('rvar_',var1,2:n_ocs,' == 1 - (AR_',var1,i,'*AR_',var1,i,' + ', 'CL_',var1,i,'*CL_',var1,i,' + 2*AR_',var1,i,'*CL_',var1,i,'*cor',i,')\n',
                                      'rvar_',var2,2:n_ocs,' == 1 - (AR_',var2,i,'*AR_',var2,i,' + ', 'CL_',var2,i,'*CL_',var2,i,' + 2*AR_',var2,i,'*CL_',var2,i,'*cor',i,')', collapse = '\n')
  }
 
  constraints = paste0(constraints, ar_var1_con,'\n', ar_var2_con,'\n', cl_var1_con,'\n', cl_var2_con, '\n', cor_con,'\n',rvar_con)
 
  f = paste0(random_intercepts, impulses, regressions, covariances, variances, constraints)
 
  if (verbose) { cat(f) }
 
  return(f)
}

[Shu Fai Cheung (張樹輝)]

I did some quick tests using the generated model you posted, using randomly generated data:

set.seed(8749314)
vars <- c("dep1", "dep2", "dep3", "dep4", "dep5", "dep6",
          "cmr1", "cmr2", "cmr3", "cmr4", "cmr5", "cmr6")
p <- length(vars)
sigma <- matrix(.2, p, p)
diag(sigma) <- 1
d <- MASS::mvrnorm(500, rep(0, p), )
colnames(d) <- vars
d <- as.data.frame(d)
head(d)

I do not yet know the source of the problem but I guess it is due to the model, not the ability of lavaan to handle noninteger exponents. I can fit Terrence's model using the CRAN version of lavaan (0.6.16). I call sem() instead of lavaan() because there may be some preset values that we need sem() to handle. Calling lavaan() is error prone because we may forget to set those values, resulting in a model different from what we expect. There are rarely cases in which there are things that can be done only by calling lavaan() but not by calling sem() (except for refitting a models using internal slots, as in some packages).

In the example above, if I run the following, with mod being the full model in the first post: 

m <- lavaan::lavaan( mod,
                      data = dat,

                      se = 'robust',
                      missing = 'fiml', # group = group,
                      meanstructure = TRUE,
                      int.ov.free = TRUE)

I got this error:

Error in sigma.inv * TT : non-conformable arrays

Not sure if this can help but you can try to add verbose = TRUE to see where the estimation failed. If I run the following using the random data:

m0 <- lavaan::sem(mod,
                  data = dat,
                  meanstructure = TRUE,
                  verbose = TRUE)

The failure occurred at attempt 4:

attempt 4 -- optim.parscale = "standardized" + start = "simple"

In your case using your real data, the output may be different. In any case, the printout may give some hints on the source of the failure.in estimation.

You can also fit a simpler model with only three waves, and see whether the problem still occurs. If estimation still fails quickly, it may be easier to diagnose the problem using a three-wave model. 

Hope this helps.

-- Shu Fai