fixing error variances from 0 to 1
et
Sunthud Pornprasertmanit
--
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 http://groups.google.com/group/lavaan.
For more options, visit https://groups.google.com/d/optout.
et
Sunthud Pornprasertmanit
model1 <- '
f1 =~ x1 + x2 + x3
x1 ~~ var1*x1
x2 ~~ var2*x2
x3 ~~ var3*x3
var1 > 0
var1 < 1
var2 > 0
var2 < 0.5
var3 > 0
var3 < 1.5
'
fit <- cfa(model1, HolzingerSwineford1939)
et
Lewis Anderson
Dear all,
Apologies for posting in an old thread - my query relates very directly to the posts here.
I am also trying to constrain an indicator's error variance to >0. However using code analogous to that posted above does not seem to be working.
(I have tried constraining only the variance of the indicator with negative variance, as shown below, I have also tried putting the same constraint on all indicators' variances. If the code above implies different constraints should be set for different indicators - e.g.
var2 > 0
var2 < 0.5
var3 > 0
var3 < 1.5
- then I do not understand why, but would like to learn.)
Perhaps I am being naïve, but I am fairly confident the model is not misspecified. To explain: the factors here are dimensions of adolescent mental wellbeing (conduct problems, emotional symptoms, hyperactivity), and the items here are teacher-rated categorical indicators (0 = 'Does not apply', 1 = 'Applies somewhat', 2 = 'Certainly applies'). In the case in which I fit the same model, except that the items are parent-rated (otherwise the same), there is no Heywood case such as occurs here. N= just below 12,000 in each case, the samples mostly overlap, there are about 9,000 cases which are in both samples.
Annotated R/lavaan code and output below. I show my original model, which showed a negative variance, my unsuccessful attempt to constrain this variance to be positive, and the analogous model with parent-rated items, for comparison.
Any help / suggestions much appreciated. More generally - the work of this group has already been a great help.
Regards,
Lewis
### Here is the model with no constraints added (OriginalModel) - to show the occurrence of the Heywood case
OriginalModel <- '
conduct =~ T3stolen + T3destroys + T3disobedient + T3lies + T3bullies + T3fights + T3irritable
emotional =~ T3worries + T3fearfulofnew + T3miserable + T3solitary
hyper =~ T3squirmy + T3restless + T3cantsettle'
fit1 <- cfa(OriginalModel, data = mydata, std.lv = TRUE, estimator = "WLSMV", ordered=c("T3stolen","T3destroys","T3disobedient","T3lies","T3bullies","T3fights","T3irritable","T3worries","T3fearfulofnew","T3miserable","T3solitary","T3squirmy","T3restless","T3cantsettle"))
### Warning message:
Warning
message:
In lav_object_post_check(lavobject) :
lavaan WARNING: some estimated
variances are negative
summary(fit1, fit.measures=TRUE)
### Output from OriginalModel / fit1 (Heywood case bold & underlined)
> summary(fit1, fit.measures=TRUE)
lavaan (0.5-20) converged normally after
23 iterations
Number of observations 11964
Estimator
DWLS Robust
Minimum Function Test Statistic 2377.605 3157.579
Degrees of freedom 74 74
P-value (Chi-square) 0.000 0.000
Scaling correction factor
0.758
Shift parameter
19.021
for simple second-order correction
(Mplus variant)
Model test baseline model:
Minimum Function
Test Statistic 283341.732 136317.028
Degrees of freedom 91 91
P-value
0.000 0.000
User model versus baseline model:
Comparative Fit
Index (CFI) 0.992 0.977
Tucker-Lewis Index (TLI) 0.990 0.972
Root Mean Square Error of Approximation:
RMSEA
0.051 0.059
90 Percent Confidence Interval 0.049
0.053 0.057 0.061
P-value RMSEA <= 0.05 0.170 0.000
Weighted Root Mean Square Residual:
WRMR 4.470 4.470
Parameter Estimates:
Information
Expected
Standard Errors Robust.sem
Latent
Variables:
Estimate Std.Err
Z-value P(>|z|)
conduct =~
T3stolen 0.786 0.011
68.537 0.000
T3destroys 0.888
0.007 129.203 0.000
T3disobedient 0.924
0.004 225.605 0.000
T3lies 0.906 0.005
182.072 0.000
T3bullies 0.899 0.006
149.104 0.000
T3fights 0.917 0.005
180.252 0.000
T3irritable 0.826
0.007 125.140 0.000
emotional =~
T3worries 0.609 0.011
55.638 0.000
T3fearfulofnew 0.723
0.012 60.289 0.000
T3miserable 1.011
0.012 84.866 0.000
T3solitary 0.580
0.013 45.556 0.000
hyper =~
T3squirmy 0.908 0.005
186.862 0.000
T3restless 0.936
0.004 224.915 0.000
T3cantsettle 0.896
0.006 158.355 0.000
Covariances:
Estimate Std.Err
Z-value P(>|z|)
conduct ~~
emotional 0.449 0.013
34.838 0.000
hyper 0.831 0.007
121.423 0.000
emotional ~~
hyper 0.452 0.013
34.680 0.000
Intercepts:
Estimate Std.Err
Z-value P(>|z|)
T3stolen 0.000
T3destroys 0.000
T3disobedient 0.000
T3lies 0.000
T3bullies 0.000
T3fights 0.000
T3irritable 0.000
T3worries 0.000
T3fearfulofnew 0.000
T3miserable 0.000
T3solitary 0.000
T3squirmy 0.000
T3restless 0.000
T3cantsettle 0.000
conduct 0.000
emotional 0.000
hyper 0.000
Thresholds:
Estimate Std.Err
Z-value P(>|z|)
T3stolen|t1 1.652
0.019 85.100 0.000
T3stolen|t2 1.996
0.025 79.326 0.000
T3destroys|t1 1.544
0.018 85.277 0.000
T3destroys|t2 2.225
0.031 72.001 0.000
T3disobednt|t1 0.957
0.014 70.447 0.000
T3disobednt|t2 1.762
0.021 84.060 0.000
T3lies|t1 1.034 0.014
73.925 0.000
T3lies|t2 1.843 0.022
82.790 0.000
T3bullies|t1 1.458
0.017 84.793 0.000
T3bullies|t2 2.083
0.027 76.811 0.000
T3fights|t1 1.348
0.016 83.349 0.000
T3fights|t2 2.117
0.028 75.729 0.000
T3irritable|t1 0.853
0.013 65.077 0.000
T3irritable|t2 1.647
0.019 85.122 0.000
T3worries|t1 0.343
0.012 29.270 0.000
T3worries|t2 1.647
0.019 85.126 0.000
T3fearflfnw|t1 0.755
0.013 59.323 0.000
T3fearflfnw|t2 1.963
0.024 80.174 0.000
T3miserable|t1 1.073
0.014 75.517 0.000
T3miserable|t2 2.054
0.026 77.673 0.000
T3solitary|t1 0.640
0.012 51.817 0.000
T3solitary|t2 1.741
0.021 84.326 0.000
T3squirmy|t1 0.985
0.014 71.756 0.000
T3squirmy|t2 1.903
0.023 81.576 0.000
T3restless|t1 0.871
0.013 66.075 0.000
T3restless|t2 1.776
0.021 83.870 0.000
T3cantsettl|t1 0.848
0.013 64.793 0.000
T3cantsettl|t2 1.771
0.021 83.940 0.000
Variances:
Estimate Std.Err
Z-value P(>|z|)
T3stolen 0.382
T3destroys 0.211
T3disobedient 0.147
T3lies 0.178
T3bullies 0.191
T3fights 0.159
T3irritable 0.318
T3worries 0.630
T3fearfulofnew 0.478
T3miserable
-0.022
T3solitary 0.664
T3squirmy 0.175
T3restless 0.124
T3cantsettle 0.196
conduct 1.000
emotional 1.000
hyper 1.000
Scales
y*:
Estimate Std.Err
Z-value P(>|z|)
T3stolen 1.000
T3destroys 1.000
T3disobedient 1.000
T3lies 1.000
T3bullies 1.000
T3fights 1.000
T3irritable 1.000
T3worries 1.000
T3fearfulofnew 1.000
T3miserable 1.000
T3solitary 1.000
T3squirmy 1.000
T3restless 1.000
T3cantsettle 1.000
### And here is my attempt to constrain the variance of T3miserable to be positive
ConstrainedModel <- '
conduct =~ T3stolen + T3destroys + T3disobedient + T3lies + T3bullies + T3fights + T3irritable
emotional =~ T3worries + T3fearfulofnew + T3miserable + T3solitary
hyper =~ T3squirmy + T3restless + T3cantsettle
T3miserable ~~ var1*T3miserable
var1 > 0
var1 < 1
'
fit2 <- cfa(ConstrainedModel, data = mydata, std.lv = TRUE, estimator = "WLSMV", ordered=c("T3stolen","T3destroys","T3disobedient","T3lies","T3bullies","T3fights","T3irritable","T3worries","T3fearfulofnew","T3miserable","T3solitary","T3squirmy","T3restless","T3cantsettle"))
### Warning messages
Warning
messages:
1: In lav_partable_constraints_ciq(partable, con = LIST, debug = debug)
:
lavaan WARNING: non-free parameter(s)
in inequality constraint(s): var1
2: In lav_object_post_check(lavobject) :
lavaan WARNING: some estimated
variances are negative
summary(fit2, fit.measures=TRUE)
### Output from ConstrainedModel / fit2 (Heywood case bold & underlined)
|
> summary(fit2, fit.measures=TRUE) Number of observations 11964 Estimator
DWLS Robust Model test baseline model: Minimum
Function Test Statistic
283341.732 136317.028 User model versus baseline model: Comparative
Fit Index (CFI)
0.992 0.977 Root Mean Square Error of Approximation: RMSEA
0.051 0.059 Weighted Root Mean Square Residual: WRMR 4.470 4.470 Parameter Estimates: Information
Expected Latent
Variables: Covariances: Intercepts: Thresholds: Variances: Scales
y*: Constraints: |
|
> |
### For comparison, here is the same again with parent-rated items
ParentModel <- '
conduct =~ P3steals + P3destroys + P3disobedient + P3lies + P3bullies + P3fights + P3irritable
emotional =~ P3worries + P3fearfulofnew + P3miserable + P3solitary
hyper =~ P3squirmy + P3restless + P3cantsettle'
fit3 <- cfa(ParentModel, data = myparentdata, std.lv = TRUE, estimator = "WLSMV", ordered=c("P3steals","P3destroys","P3disobedient","P3lies","P3bullies","P3fights","P3irritable","P3worries","P3fearfulofnew","P3miserable","P3solitary","P3squirmy","P3restless","P3cantsettle"))
### No warning messages; output below
summary(fit3, fit.measures=TRUE)
lavaan (0.5-20) converged normally after 13 iterations
Number of observations 11443
Estimator
DWLS Robust
Minimum Function Test Statistic 1954.584 2421.834
Degrees of freedom 74 74
P-value (Chi-square) 0.000 0.000
Scaling correction factor
0.810
Shift parameter
9.775
for simple second-order correction
(Mplus variant)
Model test baseline model:
Minimum Function
Test Statistic 52314.003 34660.952
Degrees of freedom 91 91
P-value
0.000 0.000
User model versus baseline model:
Comparative Fit
Index (CFI) 0.964 0.932
Tucker-Lewis Index (TLI) 0.956 0.916
Root Mean Square Error of Approximation:
RMSEA
0.047 0.053
90 Percent Confidence Interval 0.045
0.049 0.051 0.054
P-value RMSEA <= 0.05 0.995 0.007
Weighted Root Mean Square Residual:
WRMR 4.053 4.053
Parameter Estimates:
Information
Expected
Standard Errors Robust.sem
Latent
Variables:
Estimate Std.Err
Z-value P(>|z|)
conduct =~
P3steals 0.664 0.019
34.239 0.000
P3destroys 0.751
0.020 37.767 0.000
P3disobedient 0.758
0.009 82.347 0.000
P3lies 0.769 0.010
74.504 0.000
P3bullies 0.695 0.015
44.989 0.000
P3fights 0.723 0.011
65.173 0.000
P3irritable 0.704
0.009 74.766 0.000
emotional =~
P3worries 0.599 0.013
45.747 0.000
P3fearfulofnew 0.523
0.015 35.340 0.000
P3miserable 0.906
0.016 57.046 0.000
P3solitary 0.296
0.015 19.984 0.000
hyper =~
P3squirmy 0.873 0.009
92.230 0.000
P3restless 0.842
0.009 94.750 0.000
P3cantsettle 0.799
0.011 69.937 0.000
Covariances:
Estimate Std.Err
Z-value P(>|z|)
conduct ~~
emotional 0.510 0.014
35.870 0.000
hyper 0.575 0.013
43.355 0.000
emotional ~~
hyper 0.408 0.016
25.706 0.000
Intercepts:
Estimate Std.Err
Z-value P(>|z|)
P3steals 0.000
P3destroys 0.000
P3disobedient 0.000
P3lies 0.000
P3bullies 0.000
P3fights 0.000
P3irritable 0.000
P3worries 0.000
P3fearfulofnew 0.000
P3miserable 0.000
P3solitary 0.000
P3squirmy 0.000
P3restless 0.000
P3cantsettle 0.000
conduct 0.000
emotional 0.000
hyper 0.000
Thresholds:
Estimate Std.Err
Z-value P(>|z|)
P3steals|t1 1.829
0.023 81.215 0.000
P3steals|t2 2.741
0.055 49.441 0.000
P3destroys|t1 1.948
0.025 78.765 0.000
P3destroys|t2 2.548
0.044 57.662 0.000
P3disobednt|t1 0.849
0.013 63.414 0.000
P3disobednt|t2 1.965
0.025 78.370 0.000
P3lies|t1 1.132 0.015
75.968 0.000
P3lies|t2 2.198 0.031
71.374 0.000
P3bullies|t1 1.569
0.019 83.434 0.000
P3bullies|t2 2.392
0.037 64.100 0.000
P3fights|t1 1.159
0.015 76.885 0.000
P3fights|t2 2.130
0.029 73.657 0.000
P3irritable|t1 0.267
0.012 22.456 0.000
P3irritable|t2 1.235
0.016 79.063 0.000
P3worries|t1 0.279
0.012 23.499 0.000
P3worries|t2 1.415
0.017 82.489 0.000
P3fearflfnw|t1 0.776
0.013 59.293 0.000
P3fearflfnw|t2 1.768
0.022 82.133 0.000
P3miserable|t1 1.006
0.014 71.088 0.000
P3miserable|t2 1.999
0.026 77.494 0.000
P3solitary|t1 0.230
0.012 19.437 0.000
P3solitary|t2 1.105
0.015 75.023 0.000
P3squirmy|t1 1.118
0.015 75.485 0.000
P3squirmy|t2 1.698
0.020 82.897 0.000
P3restless|t1 0.729
0.013 56.439 0.000
P3restless|t2 1.426
0.017 82.607 0.000
P3cantsettl|t1 1.167
0.015 77.116 0.000
P3cantsettl|t2 1.811
0.022 81.498 0.000
Variances:
Estimate Std.Err
Z-value P(>|z|)
P3steals 0.559
P3destroys 0.436
P3disobedient 0.425
P3lies 0.409
P3bullies 0.517
P3fights 0.477
P3irritable 0.505
P3worries 0.641
P3fearfulofnew 0.726
P3miserable 0.179
P3solitary 0.912
P3squirmy 0.239
P3restless 0.291
P3cantsettle 0.362
conduct 1.000
emotional 1.000
hyper 1.000
Scales
y*:
Estimate Std.Err
Z-value P(>|z|)
P3steals 1.000
P3destroys 1.000
P3disobedient 1.000
P3lies 1.000
P3bullies 1.000
P3fights 1.000
P3irritable 1.000
P3worries 1.000
P3fearfulofnew 1.000
P3miserable 1.000
P3solitary 1.000
P3squirmy 1.000
P3restless 1.000
P3cantsettle 1.000
Lewis Anderson
Yves Rosseel
residual variances of categorical indicators. You need the 'theta'
parameterization for this.
Look for posts about theta/delta parameterization in the archive to
learn more.
Yves.
> http://groups.google.com/__group/lavaan
> <http://groups.google.com/group/lavaan>.
> For more options, visit
> https://groups.google.com/d/__optout
> <https://groups.google.com/d/optout>.
>
> --
> You received this message because you are subscribed to the
> Google Groups "lavaan" group.
> To unsubscribe from this group and stop receiving emails from
>
> --
> 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
> Visit this group at https://groups.google.com/group/lavaan.
alodie
Dear all,
in order to avoid heywood cases, I also would like to fix the error variances of the observed variables so that they are positive (i.e., > 0).
First, I tried the first solution proposed above with var > 0 and var < 1 with the following model syntax:
m1 <- 'inhib =~ antisaccade + stopsignal + colorstroop + simon + numberstroop + arrowflanker
# constraints
antisaccade ~~ b1*antisaccade
stopsignal ~~ b2*stopsignal
colorstroop ~~ b3*colorstroop
simon ~~ b4*simon
numberstroop ~~ b5*numberstroop
arrowflanker ~~ b6*arrowflanker
b1 > 0
b2 > 0
b3 > 0
b4 > 0
b5 > 0
b6 > 0
b1 < 1
b2 < 1
b3 < 1
b4 < 1
b5 < 1
b6 < 1'
fit <- sem(m1, data=d, std.lv = TRUE)
This model had bad fit statistics and strange indices (e.g., CFI = 0, RMSEA = 9.92, SRMR = 0.47), and all estimates for the error variances were 1.
Second, I tried a model with only var > 0, with the following code:
m2 <- 'inhib =~ antisaccade + stopsignal + colorstroop + simon + numberstroop + arrowflanker
# constraints
antisaccade ~~ b1*antisaccade
stopsignal ~~ b2*stopsignal
colorstroop ~~ b3*colorstroop
simon ~~ b4*simon
numberstroop ~~ b5*numberstroop
arrowflanker ~~ b6*arrowflanker
b1 > 0
b2 > 0
b3 > 0
b4 > 0
b5 > 0
b6 > 0'
fit <- sem(m2, data=d, std.lv = TRUE)
This model had good fit statistics, but with the warning that some estimated ov variances are negative. Therefore, it doesn’t seem to work.
Third, I tried a model in which all error variances were fixed to a particular value (i.e., 0.001).
m2 <- 'inhib =~ antisaccade + stopsignal + colorstroop + simon + numberstroop + arrowflanker
# constraints
antisaccade ~~ 0.001
stopsignal ~~ 0.001
colorstroop ~~ 0.001
simon ~~ 0.001
numberstroop ~~ 0.001
arrowflanker ~~ 0.001'
fit <- sem(m3, data=d, std.lv = TRUE)
But, logically, I got the error message that they were missing observed variables in dataset: ..constant..
As my data are not categorical, and the parametrization is currently used only for categorical data, this doesn’t seem to be the solution either.
Therefore: Is there a way to constraint all variances of the observed variables to be positive or at least to fix them to a particular value?
Thanks in advance for your response!
Alodie
kma...@aol.com
(1) The constraints are on the unstandardized solution. So, there is no guarantee that the unique variances are less than one. So, one likely hypothesis is that they are all greater than 1 and restricting them to 1 in the first model led them all to land at the ceiling. In this case, the fit measures may not be odd at all, but rather accurately reflecting the horrible fit to the data (not noticeably better than the baseline model).
(2) Attempting to use constraints to avoid Heywood cases is a controversial strategy (it was one issue that separated LISREL from EQS back in the day). The reason is that Heywood cases often result from model misspecification. If you constrain the estimates, the misspecification does not go away. Your problematic estimates of latent variable variances in the second model may simply result from the same underlying misspecification manifesting itself in a new way.
(3) If you want to use constraints to avoid Heywood cases, the canonical EQS strategy was to fix the unique variances to unity but freely estimate the loading for the unique variance. The loading is not bounded at either end. So, everything is in bounds. Lavaan is based on Mplus which in turn was based on LISREL. So, each of these inherits the LISREL strategy of fixing the loading and estimating the variance. Therefore, to adopt the EQS strategy, you need to not only fix the unique variance (fix, not constrain) but also to introduce a loading parameter because the LISREL model makes these implicit. One way to do this is to fix the default unique variances to zero and introduce new latent variables with unit variances to serve as unique variances.
(4) When you fix a variance to a specific value, you still need to multiply that value by the variable name on the right hand side of the ~~ operator. You omitted that in your code for the third model.
HTH,
Keith
------------------------
Keith A. Markus
John Jay College of Criminal Justice, CUNY
http://jjcweb.jjay.cuny.edu/kmarkus
Frontiers of Test Validity Theory: Measurement, Causation and Meaning.
http://www.routledge.com/books/details/9781841692203/
alodie
Hi Keith,
Thank you so much for your response. I need some clarifications about one thing. You proposed to fix the default unique variances to zero and introduce new latent variables with unit variances. If I understand you correctly, it means that the syntax of the model would be:
m4 <- 'inhib =~ antisaccade + stopsignal + colorstroop + simon + numberstroop + arrowflanker
inhib2 =~ inhib
# constraints
antisaccade ~~ 0*antisaccade
stopsignal ~~ 0*stopsignal
colorstroop ~~ 0*colorstroop
simon ~~ 0*simon
numberstroop ~~ 0*numberstroop
arrowflanker ~~ 0*arrowflanker'
fit <- sem(m4, data=d, std.lv = TRUE)
As the variance of the new latent variable “inhib2” is fixed to 1 with “std.lv = TRUE” within the sem, I did not add “inhib2~~1*inhib2” in the model syntax.
Is the syntax correct? Because when I ran the model, I got the following error message: “Error in lav_model_estimate(lavmodel = lavmodel, lavsamplestats = lavsamplestats: lavaan ERROR: initial model-implied matrix (Sigma) is not positive definite; check your model and/or starting parameters.”.
Best,
Alodie
kma...@aol.com
Alodie,
When you see that error, you can use either of the following (or both) to investigate your model specification.
lavaanify(fit)
lavInspect(fit, what='free')
You missed some steps in implementing the EQS strategy. The df should match your original model because you are just moving parameters around.
# The Holzinger and Swineford (1939) example
HS.model <- '
visual =~ 1*x1 + x2 + x3 # Define common factor
x1 ~~ 0*x1 # Zero out default unique variances
x2 ~~ 0*x2
x3 ~~ 0*x3
u1 =~ 0 # Define new latent variables for unique variances
u2 =~ 0
u3 =~ 0
visual ~~ visual # factor variance freely estimated
u1 ~~ 1*u1 # unique variances fixed to one
u2 ~~ 1*u2
u3 ~~ 1*u3
x1 ~ start(1)*u1 # Estimate loadings for unique variances
x2 ~ start(1)*u2 # The need for start values may be sample specific
x3 ~ start(1)*u3
' # end model syntax
fit <- lavaan(HS.model, data=HolzingerSwineford1939)
summary(fit)
Remember that constraints will not help you find the cause of the Heywood case. Look through the list archive for various posts by Terrence Jorgensen describing ruling out the hypothesis that the Heywood case is due to sampling error.
The EQS type model is very similar to your 2nd model constraining the unique variances to be positive. The difference is that this model does the equivalent of constraining them to be non-negative. The chi-square values should be nearly identical.