blavaan bug limiting number of variables?
188 views
Skip to first unread message
gwern branwen
Apr 25, 2016, 4:48:01 PM4/25/16
to blavaan
I have run into a problem in blavaan in trying to fit a SEM of one of my datasets. After a lot of messing with it, the problem seems to be that an internal function in blavaan, 'coeffun', errors out when you try to fit a SEM with more than 18 latent+measured variables (eg 2 latent variables made up of 8 variables each or 2+2*8=18 variables total will work, but 2 of 8 & 9 will not since 2+8+9=19; 3 latent variables of 5 each or 3+3*5=18 will work, but not 3 of 5/5/6 since 3+5+5+6=19). Jags models routinely have scores or hundreds of nodes, and lavaan seems able to estimate SEMs with dozens of variables, so I think this may be a bug in blavaan.
Some background on my system:
$ uname -a
Linux craft 3.16.0-70-generic #90~14.04.1-Ubuntu SMP Wed Apr 6 22:56:34 UTC 2016 x86_64 x86_64 x86_64 GNU/Linux
$ R
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes"
...
R> sessionInfo()
R version 3.2.5 (2016-04-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS
locale:
[1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C LC_TIME=en_US.utf8 LC_COLLATE=en_US.utf8 LC_MONETARY=en_US.utf8 LC_MESSAGES=en_US.utf8 LC_PAPER=en_US.utf8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] blavaan_0.1-3 lavaan_0.5-20 runjags_2.0.3-2 colorout_1.1-1
loaded via a namespace (and not attached):
[1] matrixStats_0.50.1 mvtnorm_1.0-5 quadprog_1.5-5 lattice_0.20-33 zoo_1.7-12 MASS_7.3-44 nonnest2_0.3 grid_3.2.5 MatrixModels_0.4-1
[10] stats4_3.2.5 coda_0.18-1 SparseM_1.7 rjags_4-4 Matrix_1.2-5 pbivnorm_0.6.0 modeest_2.1 sandwich_2.3-4 tools_3.2.5
[19] loo_0.1.6 CompQuadForm_1.4.1 parallel_3.2.5 mnormt_1.5-4 mcmc_0.9-4 MCMCpack_1.3-6 quantreg_5.21
I can successfully fit the sample PoliticalDemocracy model (it has 3 latent variables with a total of 11 measurement variables so it fits under the 18 variable limit).
The problem comes with my own dataset.
My dataset is a very large and comprehensive set of data about myself I have collected over the years. I want to model a postulated 'productivity' factor, environmental influences on it, and the effects of various randomized & nonrandomized interventions I have done. A copy of it: https://www.dropbox.com/s/8iawa1t6mwvw65t/2016-04-23-combined.csv
But when loaded up and my first stab at a SEM relating many of the variables, blavaan errors out:
R> combined <- read.csv("2016-04-23-combined.csv", header=TRUE, colClasses=c("Date", rep(NA, 101)))
R> library(blavaan)
Loading required package: runjags
Loading required package: lavaan
This is lavaan 0.5-20
lavaan is BETA software! Please report any bugs.
This is blavaan 0.1-3
blavaan is more BETA than lavaan!
R> modelS <- '
+ # Environment latent variables
+ ## Sleep quality
+ Sleep =~ ZQ + Total.Z + Time.to.Z + Time.in.Wake + Time.in.REM + Time.in.Light + Time.in.Deep + Awakenings + Morning.Feel
+ ## Room brightness
+ Light =~ Webcam.Brightness + Webcam.Color.L + Webcam.Color.A + Webcam.Color.B + CloudCover + DayLength
+ ## heat/rain/pressure
+ Weather =~ Max.TemperatureF + Mean.TemperatureF + Min.TemperatureF + Mean.Humidity + Mean.Sea.Level.PressureIn + Mean.Wind.SpeedMPH + PrecipitationIn
+
+ # unproductive latent variables
+ Music =~ Music.new + Time.Rec_Music_New + Time.Rec_Music_Old
+ IRC =~ IRC.line.count + Time.Rec_IRC
+ # master:
+ Rec =~ Music + Time.Rec_Video + Time.Rec_Fiction + Time.Rec + Time.Rec_Game + Media.consumption
+
+ # productive latent variables
+ Mnemosyne =~ Time.SRS + Mnemosyne.cards.log + Thinking.time.log + Mnemosyne.score
+ Gwern.net =~ Gwernnet.pageviews + Gwernnet.sessions + Gwernnet.visitors + Gwernnet.sessionDuration + Gwernnet.avgSessionDuration + Gwernnet.bounceRate
+ Writing =~ Patch.count.log + Gwern.net + Time.Statistics + Time.Rec_Nonfiction + Time.Writing + Time.Programming
+ Email =~ Time.Mail + Emails.sent + Emails.received
+ Bitcoin =~ Time.Bitcoin + Bitcoin.price.log + Bitcoin.30DayVolatility
+ Fitness =~ Workout.previous + Weight + Weight.BMI + Weight.body.fat + Weight.muscle + Weight.resting.metabolism + Weight.visceral.fat
+ '
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Finished running the simulation
Error in rjob$summary$statistics :
$ operator is invalid for atomic vectors
Calls: bsem -> eval -> eval -> blavaan -> coeffun
In addition: Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
R> traceback()
No traceback available
This is a reduced version of the SEM without the regressions, just the latent variables. lavaan can fit it.
Estimation output from lavaan:
R> s3 <- sem(model = modelS, se="bootstrap", bootstrap=50, missing="FIML", data = scale(combined[-1])); summary(s3)
Warning messages:
1: In lav_object_post_check(lavobject) :
lavaan WARNING: some estimated variances are negative
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
2: In lav_object_post_check(lavobject) :
lavaan WARNING: covariance matrix of latent variables is not positive definite; use inspect(fit,"cov.lv") to investigate.
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
3: In lav_object_post_check(lavobject) :
lavaan WARNING: observed variable error term matrix (theta) is not positive definite; use inspect(fit,"theta") to investigate.
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
lavaan (0.5-20) converged normally after 4064 iterations
Number of observations 2954
Number of missing patterns 65
Estimator ML
Minimum Function Test Statistic 49283.008
Degrees of freedom 1663
P-value (Chi-square) 0.000
Parameter Estimates:
Information Observed
Standard Errors Bootstrap
Number of requested bootstrap draws 50
Number of successful bootstrap draws 50
Latent Variables:
Estimate Std.Err Z-value P(>|z|)
Sleep =~
ZQ 1.000
Total.Z 1.746 0.096 18.208 0.000
Time.to.Z 0.002 0.001 1.103 0.270
Time.in.Wake 0.161 0.014 11.800 0.000
Time.in.REM 1.169 0.049 23.698 0.000
Time.in.Light 1.286 0.054 24.013 0.000
Time.in.Deep 0.666 0.033 19.888 0.000
Awakenings 0.185 0.010 18.246 0.000
Morning.Feel -0.000 0.001 -0.080 0.936
Light =~
Webcm.Brghtnss 1.000
Webcam.Color.L 1.452 0.209 6.933 0.000
Webcam.Color.A 0.193 0.032 6.073 0.000
Webcam.Color.B 0.322 0.051 6.293 0.000
CloudCover -0.107 0.031 -3.419 0.001
DayLength -0.031 0.039 -0.811 0.417
Weather =~
Max.TemperatrF 1.000
Mean.TempertrF 1.046 0.003 306.817 0.000
Min.TemperatrF 0.999 0.005 188.903 0.000
Mean.Humidity 0.150 0.017 8.936 0.000
Mn.S.Lvl.PrssI -0.197 0.013 -15.001 0.000
Men.Wnd.SpdMPH -0.117 0.019 -6.092 0.000
PrecipitatinIn -0.010 0.014 -0.705 0.481
Music =~
Music.new 1.000
Time.Rc_Msc_Nw -441.888 110.077 -4.014 0.000
Tim.Rc_Msc_Old 0.113 0.077 1.464 0.143
IRC =~
IRC.line.count 1.000
Time.Rec_IRC 0.730 0.115 6.353 0.000
Rec =~
Music 1.000
Time.Rec_Video 0.000 0.003 0.091 0.927
Time.Rec_Fictn 0.003 0.011 0.284 0.776
Time.Rec 0.008 0.027 0.309 0.758
Time.Rec_Game 0.001 0.002 0.263 0.793
Media.consmptn 0.002 0.007 0.233 0.815
Mnemosyne =~
Time.SRS 1.000
Mnmsyn.crds.lg -175.104 1161.559 -0.151 0.880
Thinking.tm.lg -267.023 1778.250 -0.150 0.881
Mnemosyne.scor 24.504 257.277 0.095 0.924
Gwern.net =~
Gwernnet.pgvws 1.000
Gwernnet.sssns 1.007 0.013 78.012 0.000
Gwernnet.vstrs 1.005 0.016 63.703 0.000
Gwrnnt.sssnDrt 0.922 0.047 19.506 0.000
Gwrnnt.vgSssnD 0.161 0.025 6.475 0.000
Gwernnet.bncRt 0.221 0.033 6.668 0.000
Writing =~
Patch.count.lg 1.000
Gwern.net 2.262 0.425 5.321 0.000
Time.Statistcs 2.525 0.648 3.895 0.000
Time.Rc_Nnfctn 2.895 0.739 3.916 0.000
Time.Writing 1.619 0.396 4.084 0.000
Time.Progrmmng 0.386 0.131 2.940 0.003
Email =~
Time.Mail 1.000
Emails.sent 0.969 0.071 13.671 0.000
Emails.receivd 1.457 0.144 10.099 0.000
Bitcoin =~
Time.Bitcoin 1.000
Bitcoin.prc.lg 173.634 72.983 2.379 0.017
Btcn.30DyVltlt 2.254 0.829 2.718 0.007
Fitness =~
Workout.previs 1.000
Weight -3.618 0.270 -13.404 0.000
Weight.BMI -3.915 0.290 -13.505 0.000
Weight.body.ft -2.321 0.192 -12.062 0.000
Weight.muscle 2.182 0.188 11.615 0.000
Wght.rstng.mtb -3.898 0.288 -13.525 0.000
Weght.vscrl.ft -3.509 0.250 -14.034 0.000
Covariances:
Estimate Std.Err Z-value P(>|z|)
Sleep ~~
Light 0.000 0.000 0.591 0.554
Weather 0.000 0.000 1.188 0.235
IRC -0.000 0.000 -1.026 0.305
Rec -0.000 0.000 -0.640 0.522
Mnemosyne -0.000 0.000 -0.031 0.975
Writing 0.000 0.000 0.073 0.942
Email -0.000 0.000 -1.163 0.245
Bitcoin -0.000 0.000 -1.181 0.238
Fitness -0.000 0.000 -0.172 0.863
Light ~~
Weather -0.041 0.014 -2.935 0.003
IRC -0.024 0.017 -1.353 0.176
Rec -0.000 0.000 -1.856 0.063
Mnemosyne 0.001 0.013 0.060 0.952
Writing 0.022 0.005 4.288 0.000
Email -0.116 0.033 -3.496 0.000
Bitcoin 0.001 0.001 2.528 0.011
Fitness 0.100 0.018 5.563 0.000
Weather ~~
IRC -0.012 0.012 -1.000 0.317
Rec -0.000 0.000 -0.481 0.630
Mnemosyne -0.000 0.003 -0.072 0.942
Writing -0.018 0.005 -3.754 0.000
Email 0.030 0.009 3.374 0.001
Bitcoin -0.000 0.000 -2.164 0.031
Fitness -0.005 0.006 -0.889 0.374
IRC ~~
Rec -0.000 0.000 -1.827 0.068
Mnemosyne 0.000 0.003 0.048 0.962
Writing 0.094 0.021 4.604 0.000
Email 0.103 0.014 7.578 0.000
Bitcoin 0.000 0.000 0.864 0.388
Fitness 0.017 0.011 1.576 0.115
Rec ~~
Mnemosyne -0.000 0.000 -0.049 0.961
Writing -0.000 0.000 -1.977 0.048
Email 0.000 0.000 1.512 0.131
Bitcoin -0.000 0.000 -1.556 0.120
Fitness -0.000 0.000 -2.039 0.041
Mnemosyne ~~
Writing 0.000 0.006 0.076 0.939
Email -0.001 0.010 -0.078 0.938
Bitcoin 0.000 0.001 0.031 0.975
Fitness 0.001 0.012 0.056 0.956
Writing ~~
Email 0.025 0.011 2.279 0.023
Bitcoin 0.001 0.000 2.628 0.009
Fitness 0.009 0.003 2.759 0.006
Email ~~
Bitcoin -0.002 0.001 -2.572 0.010
Fitness -0.107 0.011 -10.018 0.000
Bitcoin ~~
Fitness 0.001 0.001 2.477 0.013
Intercepts:
Estimate Std.Err Z-value P(>|z|)
ZQ 0.000 0.024 0.020 0.984
Total.Z 0.001 0.022 0.038 0.970
Time.to.Z 0.000 0.028 0.000 1.000
Time.in.Wake 0.000 0.023 0.003 0.997
Time.in.REM 0.001 0.026 0.022 0.982
Time.in.Light 0.001 0.021 0.030 0.976
Time.in.Deep 0.000 0.025 0.013 0.990
Awakenings 0.000 0.023 0.004 0.997
Morning.Feel 0.000 0.026 0.000 1.000
Webcm.Brghtnss -0.091 0.030 -3.071 0.002
Webcam.Color.L -0.132 0.031 -4.308 0.000
Webcam.Color.A -0.018 0.025 -0.687 0.492
Webcam.Color.B -0.029 0.025 -1.173 0.241
CloudCover -0.000 0.019 -0.001 0.999
DayLength -0.000 0.020 -0.000 1.000
Max.TemperatrF -0.000 0.019 -0.001 1.000
Mean.TempertrF -0.000 0.019 -0.001 0.999
Min.TemperatrF -0.000 0.019 -0.001 1.000
Mean.Humidity -0.000 0.019 -0.000 1.000
Mn.S.Lvl.PrssI 0.000 0.019 0.000 1.000
Men.Wnd.SpdMPH 0.000 0.018 0.000 1.000
PrecipitatinIn 0.000 0.019 0.000 1.000
Music.new -0.002 0.031 -0.059 0.953
Time.Rc_Msc_Nw -0.204 0.137 -1.488 0.137
Tim.Rc_Msc_Old 0.000 0.021 0.002 0.998
IRC.line.count -0.062 0.027 -2.274 0.023
Time.Rec_IRC -0.048 0.032 -1.499 0.134
Time.Rec_Video 0.000 0.022 0.000 1.000
Time.Rec_Fictn 0.000 0.021 0.000 1.000
Time.Rec 0.000 0.028 0.000 1.000
Time.Rec_Game 0.000 0.024 0.000 1.000
Media.consmptn 0.000 0.021 0.000 1.000
Time.SRS -0.002 0.040 -0.051 0.960
Mnmsyn.crds.lg -0.000 0.017 -0.015 0.988
Thinking.tm.lg -0.047 0.022 -2.145 0.032
Mnemosyne.scor 0.004 0.025 0.175 0.861
Gwernnet.pgvws 0.000 0.019 0.002 0.998
Gwernnet.sssns 0.000 0.019 0.002 0.998
Gwernnet.vstrs 0.000 0.019 0.002 0.998
Gwrnnt.sssnDrt 0.000 0.021 0.002 0.999
Gwrnnt.vgSssnD 0.000 0.019 0.000 1.000
Gwernnet.bncRt -0.039 0.022 -1.769 0.077
Patch.count.lg 0.000 0.022 0.001 0.999
Time.Statistcs -0.234 0.051 -4.560 0.000
Time.Rc_Nnfctn -0.269 0.050 -5.344 0.000
Time.Writing -0.150 0.040 -3.754 0.000
Time.Progrmmng -0.036 0.026 -1.385 0.166
Time.Mail 0.214 0.032 6.576 0.000
Emails.sent -0.000 0.018 -0.012 0.990
Emails.receivd -0.000 0.019 -0.017 0.986
Time.Bitcoin -0.070 0.018 -3.938 0.000
Bitcoin.prc.lg 0.001 0.019 0.028 0.977
Btcn.30DyVltlt 0.000 0.017 0.000 1.000
Workout.previs 0.000 0.016 0.016 0.987
Weight 0.760 0.060 12.767 0.000
Weight.BMI 1.325 0.064 20.661 0.000
Weight.body.ft 0.785 0.059 13.360 0.000
Weight.muscle -0.738 0.058 -12.629 0.000
Wght.rstng.mtb 1.319 0.065 20.419 0.000
Weght.vscrl.ft 1.188 0.058 20.620 0.000
Sleep 0.000
Light 0.000
Weather 0.000
Music 0.000
IRC 0.000
Rec 0.000
Mnemosyne 0.000
Gwern.net 0.000
Writing 0.000
Email 0.000
Bitcoin 0.000
Fitness 0.000
Variances:
Estimate Std.Err Z-value P(>|z|)
ZQ 0.598 0.027 22.089 0.000
Total.Z -0.237 0.014 -16.488 0.000
Time.to.Z 0.999 0.036 27.725 0.000
Time.in.Wake 0.989 0.102 9.649 0.000
Time.in.REM 0.445 0.022 20.477 0.000
Time.in.Light 0.328 0.015 21.231 0.000
Time.in.Deep 0.820 0.033 24.721 0.000
Awakenings 0.986 0.030 32.863 0.000
Morning.Feel 0.999 0.041 24.512 0.000
Webcm.Brghtnss 0.340 0.072 4.741 0.000
Webcam.Color.L -0.390 0.207 -1.884 0.060
Webcam.Color.A 0.975 0.028 35.218 0.000
Webcam.Color.B 0.931 0.036 25.722 0.000
CloudCover 0.992 0.018 56.511 0.000
DayLength 0.999 0.015 66.987 0.000
Max.TemperatrF 0.058 0.001 41.023 0.000
Mean.TempertrF -0.030 0.001 -36.718 0.000
Min.TemperatrF 0.061 0.001 40.852 0.000
Mean.Humidity 0.978 0.028 35.350 0.000
Mn.S.Lvl.PrssI 0.963 0.031 31.570 0.000
Men.Wnd.SpdMPH 0.987 0.037 26.489 0.000
PrecipitatinIn 1.000 0.035 28.538 0.000
Music.new 0.996 0.023 44.145 0.000
Time.Rc_Msc_Nw 148.624 40.568 3.664 0.000
Tim.Rc_Msc_Old 0.999 0.031 32.332 0.000
IRC.line.count 0.393 0.406 0.967 0.334
Time.Rec_IRC 0.705 0.099 7.116 0.000
Time.Rec_Video 0.999 0.039 25.528 0.000
Time.Rec_Fictn 0.999 0.054 18.578 0.000
Time.Rec 0.999 0.168 5.939 0.000
Time.Rec_Game 0.999 0.034 29.754 0.000
Media.consmptn 1.000 0.171 5.848 0.000
Time.SRS 1.000 0.039 25.420 0.000
Mnmsyn.crds.lg 0.824 0.037 22.126 0.000
Thinking.tm.lg 0.651 0.091 7.162 0.000
Mnemosyne.scor 0.997 0.045 22.007 0.000
Gwernnet.pgvws 0.013 0.002 7.133 0.000
Gwernnet.sssns -0.002 0.000 -6.271 0.000
Gwernnet.vstrs 0.003 0.000 6.679 0.000
Gwrnnt.sssnDrt 0.161 0.027 5.974 0.000
Gwrnnt.vgSssnD 0.974 0.043 22.808 0.000
Gwernnet.bncRt 0.945 0.093 10.117 0.000
Patch.count.lg 0.952 0.026 36.325 0.000
Time.Statistcs 0.785 0.070 11.177 0.000
Time.Rc_Nnfctn 0.717 0.046 15.568 0.000
Time.Writing 0.911 0.047 19.407 0.000
Time.Progrmmng 0.994 0.126 7.921 0.000
Time.Mail 0.811 0.090 8.975 0.000
Emails.sent 0.792 0.075 10.634 0.000
Emails.receivd 0.532 0.022 23.799 0.000
Time.Bitcoin 1.016 0.172 5.903 0.000
Bitcoin.prc.lg -18.607 7.592 -2.451 0.014
Btcn.30DyVltlt 0.996 0.033 30.590 0.000
Workout.previs 0.875 0.040 22.130 0.000
Weight 0.004 0.001 6.754 0.000
Weight.BMI 0.002 0.000 4.745 0.000
Weight.body.ft 0.648 0.085 7.652 0.000
Weight.muscle 0.689 0.033 20.630 0.000
Wght.rstng.mtb 0.010 0.002 6.822 0.000
Weght.vscrl.ft 0.198 0.016 12.681 0.000
Sleep 0.407 0.056 7.244 0.000
Light 0.678 0.076 8.870 0.000
Weather 0.941 0.020 48.130 0.000
Music 0.079 6.809 0.012 0.991
IRC 0.606 0.196 3.094 0.002
Rec -0.080 6.809 -0.012 0.991
Mnemosyne 0.000 0.007 0.001 0.999
Gwern.net 0.744 0.146 5.092 0.000
Writing 0.047 0.016 3.021 0.003
Email 0.220 0.034 6.415 0.000
Bitcoin 0.001 0.000 1.696 0.090
Fitness 0.125 0.015 8.244 0.000
Clearly the SEM or data are imperfect, but the output are still meaningful. Some of the variables I've found, using lavaan, have serious skew which damages the fit, and lavaan works much better if they are log-transformed. This turns out to not fix the blavaan problem:
R> combined$Time.Rec_Music_New <- log1p(combined$Time.Rec_Music_New)
R> combined$Time.Rec_Music_Old <- log1p(combined$Time.Rec_Music_Old)
R> combined$Time.Rec_Video <- log1p(combined$Time.Rec_Video)
R> combined$Time.Rec_Fiction <- log1p(combined$Time.Rec_Fiction)
R> combined$Time.Rec <- log1p(combined$Time.Rec)
R> combined$Time.Rec_Game <- log1p(combined$Time.Rec_Game)
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Finished running the simulation
Error in rjob$summary$statistics :
$ operator is invalid for atomic vectors
Calls: bsem -> eval -> eval -> blavaan -> coeffun
In addition: Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
R> traceback()
5: coeffun(lavpartable, res)
4: blavaan(model = modelS, test = "none", data = scale(combined[-1]),
model.type = "bsem", n.chains = 3, int.ov.free = TRUE, int.lv.free = FALSE,
auto.fix.first = TRUE, auto.fix.single = TRUE, auto.var = TRUE,
auto.cov.lv.x = TRUE, auto.cov.y = TRUE, auto.th = TRUE,
auto.delta = TRUE)
3: eval(expr, envir, enclos)
2: eval(mc, parent.frame())
1: bsem(model = modelS, test = "none", data = scale(combined[-1]))
The reason for the 'scale' call is to standardize each column and help avoid unit issues; the '-1' drops the 'Date' column, which, although it is not used in the model, still seems to cause problems in blavaan/lavaan ('Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric / Calls: bsem -> scale -> scale.default -> colMeans') so the 'Date' column has to be dropped to get anywhere.
I should mention that the 'scale(combined[-1])' is not the problem, since if I add a dummy column to the PoliticalDemocracy dataset and I replicate the 'bsem' call exactly, it continues to be fit fine:
tmp <- cbind(Test=rep(1,75), PoliticalDemocracy)
fit5 <- bsem(model = model, test="none", data=scale(tmp[-1])); summary(fit5)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 48 variables....
Note: Unable to calculate the multivariate psrf
Finished running the simulation
Found more than one class "Model" in cache; using the first, from namespace 'lavaan'
Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
blavaan (0.1-3) results of 10000 samples after 5000 adapt+burnin iterations
Number of observations 75
Number of missing patterns 1
Statistic
Value
Parameter Estimates:
Information MCMC
Standard Errors MCMC
Latent Variables:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ind60 =~
x1 1.000
x2 1.052 0.074 0.909 1.198 1.001 dnorm(0,1e-2)
x3 0.966 0.088 0.796 1.14 1.000 dnorm(0,1e-2)
...
At this point I started trying to cut down the model to see if it worked at all. I found that I could fit a single latent variable like 'Sleep' without triggering the error, and I could even fit two latent variables, 'Sleep' and 'Light' without triggering it:
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss); traceback() # sleep+light
Compiling rjags model...
Calling the simulation using the rjags method...
Note: the model did not require adaptation
Burning in the model for 4000 iterations...
|*************************************** | 78 4A
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 46 variables....
Finished running the simulation
Found more than one class "Model" in cache; using the first, from namespace 'lavaan'
Warning messages:
1: In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
2: In lav_data_full(data = data, group = group, group.label = group.label, :
lavaan WARNING: some cases are empty and will be removed:
2951 2952 2953 2954
Calls: bsem ... blavaan -> do.call -> lavaan -> lavData -> lav_data_full
blavaan (0.1-3) results of 10000 samples after 5000 adapt+burnin iterations
Used Total
Number of observations 2950 2954
Number of missing patterns 6
Statistic
Value
Parameter Estimates:
Information MCMC
Standard Errors MCMC
Latent Variables:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
Sleep =~
ZQ 1.000
Total.Z 1.037 0.007 1.022 1.051 1.003 dnorm(0,1e-2)
Time.to.Z -0.029 0.024 -0.076 0.018 1.000 dnorm(0,1e-2)
Time.in.Wake -0.096 0.024 -0.144 -0.048 1.000 dnorm(0,1e-2)
Time.in.REM 0.756 0.018 0.722 0.791 1.001 dnorm(0,1e-2)
Time.in.Light 0.899 0.014 0.871 0.924 1.000 dnorm(0,1e-2)
Time.in.Deep 0.501 0.022 0.459 0.543 1.000 dnorm(0,1e-2)
Awakenings 0.272 0.024 0.223 0.316 1.000 dnorm(0,1e-2)
Morning.Feel 0.424 0.024 0.377 0.472 1.000 dnorm(0,1e-2)
Light =~
Webcm.Brghtnss 1.000
Webcam.Color.L 1.023 0.008 1.006 1.039 1.000 dnorm(0,1e-2)
Webcam.Color.A 0.042 0.023 -0.002 0.087 1.000 dnorm(0,1e-2)
Webcam.Color.B 0.375 0.021 0.335 0.419 1.000 dnorm(0,1e-2)
CloudCover -0.112 0.022 -0.154 -0.067 1.000 dnorm(0,1e-2)
DayLength -0.283 0.022 -0.324 -0.239 1.000 dnorm(0,1e-2)
Covariances:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
Sleep ~~
Light -0.014 0.027 -0.066 0.038 1.001 dwish(iden,3)
Intercepts:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ZQ -0.002 0.023 -0.048 0.044 1.017 dnorm(0,1e-3)
Total.Z -0.002 0.023 -0.049 0.043 1.018 dnorm(0,1e-3)
Time.to.Z 0.000 0.023 -0.045 0.046 1.000 dnorm(0,1e-3)
Time.in.Wake 0.000 0.023 -0.044 0.048 1.000 dnorm(0,1e-3)
Time.in.REM -0.001 0.023 -0.046 0.045 1.009 dnorm(0,1e-3)
Time.in.Light -0.001 0.023 -0.046 0.044 1.014 dnorm(0,1e-3)
Time.in.Deep -0.001 0.023 -0.046 0.045 1.004 dnorm(0,1e-3)
Awakenings -0.001 0.023 -0.046 0.045 1.002 dnorm(0,1e-3)
Morning.Feel -0.015 0.024 -0.062 0.032 1.003 dnorm(0,1e-3)
Webcm.Brghtnss -0.020 0.022 -0.062 0.023 1.008 dnorm(0,1e-3)
Webcam.Color.L -0.020 0.022 -0.062 0.023 1.009 dnorm(0,1e-3)
Webcam.Color.A -0.001 0.021 -0.044 0.04 1.000 dnorm(0,1e-3)
Webcam.Color.B -0.007 0.021 -0.049 0.034 1.001 dnorm(0,1e-3)
CloudCover 0.000 0.018 -0.037 0.036 1.000 dnorm(0,1e-3)
DayLength -0.000 0.018 -0.035 0.037 1.000 dnorm(0,1e-3)
Sleep 0.000
Light 0.000
Variances:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ZQ 0.077 0.003 0.072 0.083 1.000 dgamma(1,.5)
Total.Z 0.009 0.001 0.008 0.011 1.001 dgamma(1,.5)
Time.to.Z 1.000 0.033 0.936 1.064 1.000 dgamma(1,.5)
Time.in.Wake 0.992 0.033 0.929 1.057 1.000 dgamma(1,.5)
Time.in.REM 0.473 0.016 0.444 0.505 1.000 dgamma(1,.5)
Time.in.Light 0.256 0.009 0.239 0.274 1.000 dgamma(1,.5)
Time.in.Deep 0.769 0.025 0.72 0.819 1.000 dgamma(1,.5)
Awakenings 0.933 0.031 0.872 0.992 1.000 dgamma(1,.5)
Morning.Feel 0.855 0.028 0.8 0.91 1.000 dgamma(1,.5)
Webcm.Brghtnss 0.069 0.005 0.06 0.079 1.001 dgamma(1,.5)
Webcam.Color.L 0.027 0.004 0.018 0.035 1.002 dgamma(1,.5)
Webcam.Color.A 0.999 0.030 0.941 1.059 1.000 dgamma(1,.5)
Webcam.Color.B 0.870 0.026 0.819 0.922 1.000 dgamma(1,.5)
CloudCover 0.989 0.026 0.94 1.04 1.000 dgamma(1,.5)
DayLength 0.926 0.025 0.879 0.977 1.000 dgamma(1,.5)
Sleep 0.924 0.033 0.861 0.988 1.000 dwish(iden,3)
Light 0.933 0.031 0.874 0.993 1.000 dwish(iden,3)
It takes several hours to fit using blavaan, which was getting in the way of trying to identify minimal cases, so I cut down on sampling drastically ('burnin=1000, sample=1000'). Hopefully that doesn't interfere with triggering the bug, but it made it a lot faster to iterate.
Tweaking the model repeatedly I found:
- Political/Democracy builtin dataset: 3 latent variables with 11 measurement variables: worked
- Sleep=2,Light=3,Weather=3,Music=3,IRC=2, 13: worked
- Sleep=3,Light=3,Weather=3,Music=3,IRC=2, 14: error
- Sleep+Light+Weather, with Weather reduced to 2 measurements instead of 7: error
- Sleep=9+Light=6+Weather=2, 17: error
- Sleep=8+Light=6+Weather=2, 16: error
- Sleep=7+Light=6+Weather=2, 15: worked
- Sleep=7+Light=6+Weather=3, 16: error
- Sleep=6+Light=6+Weather=3, 15: worked
- Sleep=6+Light=5+Weather=4, 15: worked
- Sleep=6+Light=5+Weather=5, 16: error
- Sleep=5+Light=5+Weather=5, 15: worked
- Sleep=4+Light=4+Weather=4+Music=3, 15: error
- Sleep=9,Weather=7, 16: worked
- Sleep=9,Weather=9, 18: error
- Sleep=8,Weather=9, 17: error
- Sleep=8,Weather=8, 16: worked
So with 5 latent variables, blavaan will fit models with <=13 measurement variables. With 3, it'll fit <=15 With 2 latent variables, it'll fit <=16 measurement variables. (I don't have any latent variables with >=17 measurement variables, so I couldn't check the singular case.) So 5+13, 3+15, or 2+16. This suggests a pattern of # of latent variables + # of measurement variables <= 18.
I looked at coeffun's call in 'bsem', but aside from indicating that it's called after JAGS is done running (successfully?), I didn't see anything that helped me:
...
verbose <- lavoptions$verbose
start.time <- proc.time()[3]
x <- NULL
if (lavm...@nx.free > 0L) {
if (!trans.exists) {
jagtrans <- try(lav2jags(model = lavpartable, lavdata = lavdata,
ov.cp = ov.cp, lv.cp = lv.cp, lv.x.wish = lavoptions$auto.cov.lv.x,
dp = dp, n.chains = n.chains, jagextra = jagextra,
inits = inits), silent = TRUE)
}
if (!inherits(jagtrans, "try-error")) {
if (jagfile) {
dir.create(path = jagdir, showWarnings = FALSE)
cat(jagtrans$model, file = paste(jagdir, "/sem.jag",
sep = ""))
save(jagtrans, file = paste(jagdir, "/semjags.rda",
sep = ""))
}
lavpartable <- mergejag(lavpartable, jagtrans$coefvec)
sampparms <- jagtrans$coefvec[, 1]
if ("monitor" %in% names(jagextra)) {
sampparms <- c(sampparms, jagextra$monitor)
}
if (inits == "jags")
jagtrans$inits <- NA
rjarg <- with(jagtrans, list(model = paste(model),
monitor = sampparms, data = data, inits = inits))
rjarg <- c(rjarg, mfj, jagcontrol)
if (suppressMessages(requireNamespace("modeest")))
runjags.options(mode.continuous = TRUE)
if (jag.do.fit) {
res <- try(do.call("run.jags", rjarg))
}
else {
res <- NULL
}
if (inherits(res, "try-error")) {
if (!trans.exists) {
dir.create(path = jagdir, showWarnings = FALSE)
cat(jagtrans$model, file = paste(jagdir, "/sem.jag",
sep = ""))
save(jagtrans, file = paste(jagdir, "/semjags.rda",
sep = ""))
}
stop("blavaan ERROR: problem with jags estimation. The jags model and data have been exported.")
}
}
else {
stop("blavaan ERROR: problem with translation from lavaan to jags.")
}
parests <- coeffun(lavpartable, res)
x <- parests$x
lavpartable <- parests$lavpartable
...
Some background on my system:
$ uname -a
Linux craft 3.16.0-70-generic #90~14.04.1-Ubuntu SMP Wed Apr 6 22:56:34 UTC 2016 x86_64 x86_64 x86_64 GNU/Linux
$ R
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes"
...
R> sessionInfo()
R version 3.2.5 (2016-04-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS
locale:
[1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C LC_TIME=en_US.utf8 LC_COLLATE=en_US.utf8 LC_MONETARY=en_US.utf8 LC_MESSAGES=en_US.utf8 LC_PAPER=en_US.utf8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] blavaan_0.1-3 lavaan_0.5-20 runjags_2.0.3-2 colorout_1.1-1
loaded via a namespace (and not attached):
[1] matrixStats_0.50.1 mvtnorm_1.0-5 quadprog_1.5-5 lattice_0.20-33 zoo_1.7-12 MASS_7.3-44 nonnest2_0.3 grid_3.2.5 MatrixModels_0.4-1
[10] stats4_3.2.5 coda_0.18-1 SparseM_1.7 rjags_4-4 Matrix_1.2-5 pbivnorm_0.6.0 modeest_2.1 sandwich_2.3-4 tools_3.2.5
[19] loo_0.1.6 CompQuadForm_1.4.1 parallel_3.2.5 mnormt_1.5-4 mcmc_0.9-4 MCMCpack_1.3-6 quantreg_5.21
I can successfully fit the sample PoliticalDemocracy model (it has 3 latent variables with a total of 11 measurement variables so it fits under the 18 variable limit).
The problem comes with my own dataset.
My dataset is a very large and comprehensive set of data about myself I have collected over the years. I want to model a postulated 'productivity' factor, environmental influences on it, and the effects of various randomized & nonrandomized interventions I have done. A copy of it: https://www.dropbox.com/s/8iawa1t6mwvw65t/2016-04-23-combined.csv
But when loaded up and my first stab at a SEM relating many of the variables, blavaan errors out:
R> combined <- read.csv("2016-04-23-combined.csv", header=TRUE, colClasses=c("Date", rep(NA, 101)))
R> library(blavaan)
Loading required package: runjags
Loading required package: lavaan
This is lavaan 0.5-20
lavaan is BETA software! Please report any bugs.
This is blavaan 0.1-3
blavaan is more BETA than lavaan!
R> modelS <- '
+ # Environment latent variables
+ ## Sleep quality
+ Sleep =~ ZQ + Total.Z + Time.to.Z + Time.in.Wake + Time.in.REM + Time.in.Light + Time.in.Deep + Awakenings + Morning.Feel
+ ## Room brightness
+ Light =~ Webcam.Brightness + Webcam.Color.L + Webcam.Color.A + Webcam.Color.B + CloudCover + DayLength
+ ## heat/rain/pressure
+ Weather =~ Max.TemperatureF + Mean.TemperatureF + Min.TemperatureF + Mean.Humidity + Mean.Sea.Level.PressureIn + Mean.Wind.SpeedMPH + PrecipitationIn
+
+ # unproductive latent variables
+ Music =~ Music.new + Time.Rec_Music_New + Time.Rec_Music_Old
+ IRC =~ IRC.line.count + Time.Rec_IRC
+ # master:
+ Rec =~ Music + Time.Rec_Video + Time.Rec_Fiction + Time.Rec + Time.Rec_Game + Media.consumption
+
+ # productive latent variables
+ Mnemosyne =~ Time.SRS + Mnemosyne.cards.log + Thinking.time.log + Mnemosyne.score
+ Gwern.net =~ Gwernnet.pageviews + Gwernnet.sessions + Gwernnet.visitors + Gwernnet.sessionDuration + Gwernnet.avgSessionDuration + Gwernnet.bounceRate
+ Writing =~ Patch.count.log + Gwern.net + Time.Statistics + Time.Rec_Nonfiction + Time.Writing + Time.Programming
+ Email =~ Time.Mail + Emails.sent + Emails.received
+ Bitcoin =~ Time.Bitcoin + Bitcoin.price.log + Bitcoin.30DayVolatility
+ Fitness =~ Workout.previous + Weight + Weight.BMI + Weight.body.fat + Weight.muscle + Weight.resting.metabolism + Weight.visceral.fat
+ '
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Finished running the simulation
Error in rjob$summary$statistics :
$ operator is invalid for atomic vectors
Calls: bsem -> eval -> eval -> blavaan -> coeffun
In addition: Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
R> traceback()
No traceback available
This is a reduced version of the SEM without the regressions, just the latent variables. lavaan can fit it.
Estimation output from lavaan:
R> s3 <- sem(model = modelS, se="bootstrap", bootstrap=50, missing="FIML", data = scale(combined[-1])); summary(s3)
Warning messages:
1: In lav_object_post_check(lavobject) :
lavaan WARNING: some estimated variances are negative
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
2: In lav_object_post_check(lavobject) :
lavaan WARNING: covariance matrix of latent variables is not positive definite; use inspect(fit,"cov.lv") to investigate.
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
3: In lav_object_post_check(lavobject) :
lavaan WARNING: observed variable error term matrix (theta) is not positive definite; use inspect(fit,"theta") to investigate.
Calls: sem ... eval -> <Anonymous> -> lavInspect -> lav_object_post_check
lavaan (0.5-20) converged normally after 4064 iterations
Number of observations 2954
Number of missing patterns 65
Estimator ML
Minimum Function Test Statistic 49283.008
Degrees of freedom 1663
P-value (Chi-square) 0.000
Parameter Estimates:
Information Observed
Standard Errors Bootstrap
Number of requested bootstrap draws 50
Number of successful bootstrap draws 50
Latent Variables:
Estimate Std.Err Z-value P(>|z|)
Sleep =~
ZQ 1.000
Total.Z 1.746 0.096 18.208 0.000
Time.to.Z 0.002 0.001 1.103 0.270
Time.in.Wake 0.161 0.014 11.800 0.000
Time.in.REM 1.169 0.049 23.698 0.000
Time.in.Light 1.286 0.054 24.013 0.000
Time.in.Deep 0.666 0.033 19.888 0.000
Awakenings 0.185 0.010 18.246 0.000
Morning.Feel -0.000 0.001 -0.080 0.936
Light =~
Webcm.Brghtnss 1.000
Webcam.Color.L 1.452 0.209 6.933 0.000
Webcam.Color.A 0.193 0.032 6.073 0.000
Webcam.Color.B 0.322 0.051 6.293 0.000
CloudCover -0.107 0.031 -3.419 0.001
DayLength -0.031 0.039 -0.811 0.417
Weather =~
Max.TemperatrF 1.000
Mean.TempertrF 1.046 0.003 306.817 0.000
Min.TemperatrF 0.999 0.005 188.903 0.000
Mean.Humidity 0.150 0.017 8.936 0.000
Mn.S.Lvl.PrssI -0.197 0.013 -15.001 0.000
Men.Wnd.SpdMPH -0.117 0.019 -6.092 0.000
PrecipitatinIn -0.010 0.014 -0.705 0.481
Music =~
Music.new 1.000
Time.Rc_Msc_Nw -441.888 110.077 -4.014 0.000
Tim.Rc_Msc_Old 0.113 0.077 1.464 0.143
IRC =~
IRC.line.count 1.000
Time.Rec_IRC 0.730 0.115 6.353 0.000
Rec =~
Music 1.000
Time.Rec_Video 0.000 0.003 0.091 0.927
Time.Rec_Fictn 0.003 0.011 0.284 0.776
Time.Rec 0.008 0.027 0.309 0.758
Time.Rec_Game 0.001 0.002 0.263 0.793
Media.consmptn 0.002 0.007 0.233 0.815
Mnemosyne =~
Time.SRS 1.000
Mnmsyn.crds.lg -175.104 1161.559 -0.151 0.880
Thinking.tm.lg -267.023 1778.250 -0.150 0.881
Mnemosyne.scor 24.504 257.277 0.095 0.924
Gwern.net =~
Gwernnet.pgvws 1.000
Gwernnet.sssns 1.007 0.013 78.012 0.000
Gwernnet.vstrs 1.005 0.016 63.703 0.000
Gwrnnt.sssnDrt 0.922 0.047 19.506 0.000
Gwrnnt.vgSssnD 0.161 0.025 6.475 0.000
Gwernnet.bncRt 0.221 0.033 6.668 0.000
Writing =~
Patch.count.lg 1.000
Gwern.net 2.262 0.425 5.321 0.000
Time.Statistcs 2.525 0.648 3.895 0.000
Time.Rc_Nnfctn 2.895 0.739 3.916 0.000
Time.Writing 1.619 0.396 4.084 0.000
Time.Progrmmng 0.386 0.131 2.940 0.003
Email =~
Time.Mail 1.000
Emails.sent 0.969 0.071 13.671 0.000
Emails.receivd 1.457 0.144 10.099 0.000
Bitcoin =~
Time.Bitcoin 1.000
Bitcoin.prc.lg 173.634 72.983 2.379 0.017
Btcn.30DyVltlt 2.254 0.829 2.718 0.007
Fitness =~
Workout.previs 1.000
Weight -3.618 0.270 -13.404 0.000
Weight.BMI -3.915 0.290 -13.505 0.000
Weight.body.ft -2.321 0.192 -12.062 0.000
Weight.muscle 2.182 0.188 11.615 0.000
Wght.rstng.mtb -3.898 0.288 -13.525 0.000
Weght.vscrl.ft -3.509 0.250 -14.034 0.000
Covariances:
Estimate Std.Err Z-value P(>|z|)
Sleep ~~
Light 0.000 0.000 0.591 0.554
Weather 0.000 0.000 1.188 0.235
IRC -0.000 0.000 -1.026 0.305
Rec -0.000 0.000 -0.640 0.522
Mnemosyne -0.000 0.000 -0.031 0.975
Writing 0.000 0.000 0.073 0.942
Email -0.000 0.000 -1.163 0.245
Bitcoin -0.000 0.000 -1.181 0.238
Fitness -0.000 0.000 -0.172 0.863
Light ~~
Weather -0.041 0.014 -2.935 0.003
IRC -0.024 0.017 -1.353 0.176
Rec -0.000 0.000 -1.856 0.063
Mnemosyne 0.001 0.013 0.060 0.952
Writing 0.022 0.005 4.288 0.000
Email -0.116 0.033 -3.496 0.000
Bitcoin 0.001 0.001 2.528 0.011
Fitness 0.100 0.018 5.563 0.000
Weather ~~
IRC -0.012 0.012 -1.000 0.317
Rec -0.000 0.000 -0.481 0.630
Mnemosyne -0.000 0.003 -0.072 0.942
Writing -0.018 0.005 -3.754 0.000
Email 0.030 0.009 3.374 0.001
Bitcoin -0.000 0.000 -2.164 0.031
Fitness -0.005 0.006 -0.889 0.374
IRC ~~
Rec -0.000 0.000 -1.827 0.068
Mnemosyne 0.000 0.003 0.048 0.962
Writing 0.094 0.021 4.604 0.000
Email 0.103 0.014 7.578 0.000
Bitcoin 0.000 0.000 0.864 0.388
Fitness 0.017 0.011 1.576 0.115
Rec ~~
Mnemosyne -0.000 0.000 -0.049 0.961
Writing -0.000 0.000 -1.977 0.048
Email 0.000 0.000 1.512 0.131
Bitcoin -0.000 0.000 -1.556 0.120
Fitness -0.000 0.000 -2.039 0.041
Mnemosyne ~~
Writing 0.000 0.006 0.076 0.939
Email -0.001 0.010 -0.078 0.938
Bitcoin 0.000 0.001 0.031 0.975
Fitness 0.001 0.012 0.056 0.956
Writing ~~
Email 0.025 0.011 2.279 0.023
Bitcoin 0.001 0.000 2.628 0.009
Fitness 0.009 0.003 2.759 0.006
Email ~~
Bitcoin -0.002 0.001 -2.572 0.010
Fitness -0.107 0.011 -10.018 0.000
Bitcoin ~~
Fitness 0.001 0.001 2.477 0.013
Intercepts:
Estimate Std.Err Z-value P(>|z|)
ZQ 0.000 0.024 0.020 0.984
Total.Z 0.001 0.022 0.038 0.970
Time.to.Z 0.000 0.028 0.000 1.000
Time.in.Wake 0.000 0.023 0.003 0.997
Time.in.REM 0.001 0.026 0.022 0.982
Time.in.Light 0.001 0.021 0.030 0.976
Time.in.Deep 0.000 0.025 0.013 0.990
Awakenings 0.000 0.023 0.004 0.997
Morning.Feel 0.000 0.026 0.000 1.000
Webcm.Brghtnss -0.091 0.030 -3.071 0.002
Webcam.Color.L -0.132 0.031 -4.308 0.000
Webcam.Color.A -0.018 0.025 -0.687 0.492
Webcam.Color.B -0.029 0.025 -1.173 0.241
CloudCover -0.000 0.019 -0.001 0.999
DayLength -0.000 0.020 -0.000 1.000
Max.TemperatrF -0.000 0.019 -0.001 1.000
Mean.TempertrF -0.000 0.019 -0.001 0.999
Min.TemperatrF -0.000 0.019 -0.001 1.000
Mean.Humidity -0.000 0.019 -0.000 1.000
Mn.S.Lvl.PrssI 0.000 0.019 0.000 1.000
Men.Wnd.SpdMPH 0.000 0.018 0.000 1.000
PrecipitatinIn 0.000 0.019 0.000 1.000
Music.new -0.002 0.031 -0.059 0.953
Time.Rc_Msc_Nw -0.204 0.137 -1.488 0.137
Tim.Rc_Msc_Old 0.000 0.021 0.002 0.998
IRC.line.count -0.062 0.027 -2.274 0.023
Time.Rec_IRC -0.048 0.032 -1.499 0.134
Time.Rec_Video 0.000 0.022 0.000 1.000
Time.Rec_Fictn 0.000 0.021 0.000 1.000
Time.Rec 0.000 0.028 0.000 1.000
Time.Rec_Game 0.000 0.024 0.000 1.000
Media.consmptn 0.000 0.021 0.000 1.000
Time.SRS -0.002 0.040 -0.051 0.960
Mnmsyn.crds.lg -0.000 0.017 -0.015 0.988
Thinking.tm.lg -0.047 0.022 -2.145 0.032
Mnemosyne.scor 0.004 0.025 0.175 0.861
Gwernnet.pgvws 0.000 0.019 0.002 0.998
Gwernnet.sssns 0.000 0.019 0.002 0.998
Gwernnet.vstrs 0.000 0.019 0.002 0.998
Gwrnnt.sssnDrt 0.000 0.021 0.002 0.999
Gwrnnt.vgSssnD 0.000 0.019 0.000 1.000
Gwernnet.bncRt -0.039 0.022 -1.769 0.077
Patch.count.lg 0.000 0.022 0.001 0.999
Time.Statistcs -0.234 0.051 -4.560 0.000
Time.Rc_Nnfctn -0.269 0.050 -5.344 0.000
Time.Writing -0.150 0.040 -3.754 0.000
Time.Progrmmng -0.036 0.026 -1.385 0.166
Time.Mail 0.214 0.032 6.576 0.000
Emails.sent -0.000 0.018 -0.012 0.990
Emails.receivd -0.000 0.019 -0.017 0.986
Time.Bitcoin -0.070 0.018 -3.938 0.000
Bitcoin.prc.lg 0.001 0.019 0.028 0.977
Btcn.30DyVltlt 0.000 0.017 0.000 1.000
Workout.previs 0.000 0.016 0.016 0.987
Weight 0.760 0.060 12.767 0.000
Weight.BMI 1.325 0.064 20.661 0.000
Weight.body.ft 0.785 0.059 13.360 0.000
Weight.muscle -0.738 0.058 -12.629 0.000
Wght.rstng.mtb 1.319 0.065 20.419 0.000
Weght.vscrl.ft 1.188 0.058 20.620 0.000
Sleep 0.000
Light 0.000
Weather 0.000
Music 0.000
IRC 0.000
Rec 0.000
Mnemosyne 0.000
Gwern.net 0.000
Writing 0.000
Email 0.000
Bitcoin 0.000
Fitness 0.000
Variances:
Estimate Std.Err Z-value P(>|z|)
ZQ 0.598 0.027 22.089 0.000
Total.Z -0.237 0.014 -16.488 0.000
Time.to.Z 0.999 0.036 27.725 0.000
Time.in.Wake 0.989 0.102 9.649 0.000
Time.in.REM 0.445 0.022 20.477 0.000
Time.in.Light 0.328 0.015 21.231 0.000
Time.in.Deep 0.820 0.033 24.721 0.000
Awakenings 0.986 0.030 32.863 0.000
Morning.Feel 0.999 0.041 24.512 0.000
Webcm.Brghtnss 0.340 0.072 4.741 0.000
Webcam.Color.L -0.390 0.207 -1.884 0.060
Webcam.Color.A 0.975 0.028 35.218 0.000
Webcam.Color.B 0.931 0.036 25.722 0.000
CloudCover 0.992 0.018 56.511 0.000
DayLength 0.999 0.015 66.987 0.000
Max.TemperatrF 0.058 0.001 41.023 0.000
Mean.TempertrF -0.030 0.001 -36.718 0.000
Min.TemperatrF 0.061 0.001 40.852 0.000
Mean.Humidity 0.978 0.028 35.350 0.000
Mn.S.Lvl.PrssI 0.963 0.031 31.570 0.000
Men.Wnd.SpdMPH 0.987 0.037 26.489 0.000
PrecipitatinIn 1.000 0.035 28.538 0.000
Music.new 0.996 0.023 44.145 0.000
Time.Rc_Msc_Nw 148.624 40.568 3.664 0.000
Tim.Rc_Msc_Old 0.999 0.031 32.332 0.000
IRC.line.count 0.393 0.406 0.967 0.334
Time.Rec_IRC 0.705 0.099 7.116 0.000
Time.Rec_Video 0.999 0.039 25.528 0.000
Time.Rec_Fictn 0.999 0.054 18.578 0.000
Time.Rec 0.999 0.168 5.939 0.000
Time.Rec_Game 0.999 0.034 29.754 0.000
Media.consmptn 1.000 0.171 5.848 0.000
Time.SRS 1.000 0.039 25.420 0.000
Mnmsyn.crds.lg 0.824 0.037 22.126 0.000
Thinking.tm.lg 0.651 0.091 7.162 0.000
Mnemosyne.scor 0.997 0.045 22.007 0.000
Gwernnet.pgvws 0.013 0.002 7.133 0.000
Gwernnet.sssns -0.002 0.000 -6.271 0.000
Gwernnet.vstrs 0.003 0.000 6.679 0.000
Gwrnnt.sssnDrt 0.161 0.027 5.974 0.000
Gwrnnt.vgSssnD 0.974 0.043 22.808 0.000
Gwernnet.bncRt 0.945 0.093 10.117 0.000
Patch.count.lg 0.952 0.026 36.325 0.000
Time.Statistcs 0.785 0.070 11.177 0.000
Time.Rc_Nnfctn 0.717 0.046 15.568 0.000
Time.Writing 0.911 0.047 19.407 0.000
Time.Progrmmng 0.994 0.126 7.921 0.000
Time.Mail 0.811 0.090 8.975 0.000
Emails.sent 0.792 0.075 10.634 0.000
Emails.receivd 0.532 0.022 23.799 0.000
Time.Bitcoin 1.016 0.172 5.903 0.000
Bitcoin.prc.lg -18.607 7.592 -2.451 0.014
Btcn.30DyVltlt 0.996 0.033 30.590 0.000
Workout.previs 0.875 0.040 22.130 0.000
Weight 0.004 0.001 6.754 0.000
Weight.BMI 0.002 0.000 4.745 0.000
Weight.body.ft 0.648 0.085 7.652 0.000
Weight.muscle 0.689 0.033 20.630 0.000
Wght.rstng.mtb 0.010 0.002 6.822 0.000
Weght.vscrl.ft 0.198 0.016 12.681 0.000
Sleep 0.407 0.056 7.244 0.000
Light 0.678 0.076 8.870 0.000
Weather 0.941 0.020 48.130 0.000
Music 0.079 6.809 0.012 0.991
IRC 0.606 0.196 3.094 0.002
Rec -0.080 6.809 -0.012 0.991
Mnemosyne 0.000 0.007 0.001 0.999
Gwern.net 0.744 0.146 5.092 0.000
Writing 0.047 0.016 3.021 0.003
Email 0.220 0.034 6.415 0.000
Bitcoin 0.001 0.000 1.696 0.090
Fitness 0.125 0.015 8.244 0.000
Clearly the SEM or data are imperfect, but the output are still meaningful. Some of the variables I've found, using lavaan, have serious skew which damages the fit, and lavaan works much better if they are log-transformed. This turns out to not fix the blavaan problem:
R> combined$Time.Rec_Music_New <- log1p(combined$Time.Rec_Music_New)
R> combined$Time.Rec_Music_Old <- log1p(combined$Time.Rec_Music_Old)
R> combined$Time.Rec_Video <- log1p(combined$Time.Rec_Video)
R> combined$Time.Rec_Fiction <- log1p(combined$Time.Rec_Fiction)
R> combined$Time.Rec <- log1p(combined$Time.Rec)
R> combined$Time.Rec_Game <- log1p(combined$Time.Rec_Game)
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Finished running the simulation
Error in rjob$summary$statistics :
$ operator is invalid for atomic vectors
Calls: bsem -> eval -> eval -> blavaan -> coeffun
In addition: Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
R> traceback()
5: coeffun(lavpartable, res)
4: blavaan(model = modelS, test = "none", data = scale(combined[-1]),
model.type = "bsem", n.chains = 3, int.ov.free = TRUE, int.lv.free = FALSE,
auto.fix.first = TRUE, auto.fix.single = TRUE, auto.var = TRUE,
auto.cov.lv.x = TRUE, auto.cov.y = TRUE, auto.th = TRUE,
auto.delta = TRUE)
3: eval(expr, envir, enclos)
2: eval(mc, parent.frame())
1: bsem(model = modelS, test = "none", data = scale(combined[-1]))
The reason for the 'scale' call is to standardize each column and help avoid unit issues; the '-1' drops the 'Date' column, which, although it is not used in the model, still seems to cause problems in blavaan/lavaan ('Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric / Calls: bsem -> scale -> scale.default -> colMeans') so the 'Date' column has to be dropped to get anywhere.
I should mention that the 'scale(combined[-1])' is not the problem, since if I add a dummy column to the PoliticalDemocracy dataset and I replicate the 'bsem' call exactly, it continues to be fit fine:
tmp <- cbind(Test=rep(1,75), PoliticalDemocracy)
fit5 <- bsem(model = model, test="none", data=scale(tmp[-1])); summary(fit5)
Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 1000 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 4000 iterations...
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 48 variables....
Note: Unable to calculate the multivariate psrf
Finished running the simulation
Found more than one class "Model" in cache; using the first, from namespace 'lavaan'
Warning message:
In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
blavaan (0.1-3) results of 10000 samples after 5000 adapt+burnin iterations
Number of observations 75
Number of missing patterns 1
Statistic
Value
Parameter Estimates:
Information MCMC
Standard Errors MCMC
Latent Variables:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ind60 =~
x1 1.000
x2 1.052 0.074 0.909 1.198 1.001 dnorm(0,1e-2)
x3 0.966 0.088 0.796 1.14 1.000 dnorm(0,1e-2)
...
At this point I started trying to cut down the model to see if it worked at all. I found that I could fit a single latent variable like 'Sleep' without triggering the error, and I could even fit two latent variables, 'Sleep' and 'Light' without triggering it:
R> ss <- bsem(model = modelS, test="none", data = scale(combined[-1])); summary(ss); traceback() # sleep+light
Compiling rjags model...
Calling the simulation using the rjags method...
Note: the model did not require adaptation
Burning in the model for 4000 iterations...
|*************************************** | 78 4A
|**************************************************| 100%
Running the model for 10000 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 46 variables....
Finished running the simulation
Found more than one class "Model" in cache; using the first, from namespace 'lavaan'
Warning messages:
1: In lavData(data = data, group = group, group.label = group.label, :
lavaan WARNING: data argument has been coerced to a data.frame
Calls: bsem ... eval -> eval -> blavaan -> do.call -> lavaan -> lavData
2: In lav_data_full(data = data, group = group, group.label = group.label, :
lavaan WARNING: some cases are empty and will be removed:
2951 2952 2953 2954
Calls: bsem ... blavaan -> do.call -> lavaan -> lavData -> lav_data_full
blavaan (0.1-3) results of 10000 samples after 5000 adapt+burnin iterations
Used Total
Number of observations 2950 2954
Number of missing patterns 6
Statistic
Value
Parameter Estimates:
Information MCMC
Standard Errors MCMC
Latent Variables:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
Sleep =~
ZQ 1.000
Total.Z 1.037 0.007 1.022 1.051 1.003 dnorm(0,1e-2)
Time.to.Z -0.029 0.024 -0.076 0.018 1.000 dnorm(0,1e-2)
Time.in.Wake -0.096 0.024 -0.144 -0.048 1.000 dnorm(0,1e-2)
Time.in.REM 0.756 0.018 0.722 0.791 1.001 dnorm(0,1e-2)
Time.in.Light 0.899 0.014 0.871 0.924 1.000 dnorm(0,1e-2)
Time.in.Deep 0.501 0.022 0.459 0.543 1.000 dnorm(0,1e-2)
Awakenings 0.272 0.024 0.223 0.316 1.000 dnorm(0,1e-2)
Morning.Feel 0.424 0.024 0.377 0.472 1.000 dnorm(0,1e-2)
Light =~
Webcm.Brghtnss 1.000
Webcam.Color.L 1.023 0.008 1.006 1.039 1.000 dnorm(0,1e-2)
Webcam.Color.A 0.042 0.023 -0.002 0.087 1.000 dnorm(0,1e-2)
Webcam.Color.B 0.375 0.021 0.335 0.419 1.000 dnorm(0,1e-2)
CloudCover -0.112 0.022 -0.154 -0.067 1.000 dnorm(0,1e-2)
DayLength -0.283 0.022 -0.324 -0.239 1.000 dnorm(0,1e-2)
Covariances:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
Sleep ~~
Light -0.014 0.027 -0.066 0.038 1.001 dwish(iden,3)
Intercepts:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ZQ -0.002 0.023 -0.048 0.044 1.017 dnorm(0,1e-3)
Total.Z -0.002 0.023 -0.049 0.043 1.018 dnorm(0,1e-3)
Time.to.Z 0.000 0.023 -0.045 0.046 1.000 dnorm(0,1e-3)
Time.in.Wake 0.000 0.023 -0.044 0.048 1.000 dnorm(0,1e-3)
Time.in.REM -0.001 0.023 -0.046 0.045 1.009 dnorm(0,1e-3)
Time.in.Light -0.001 0.023 -0.046 0.044 1.014 dnorm(0,1e-3)
Time.in.Deep -0.001 0.023 -0.046 0.045 1.004 dnorm(0,1e-3)
Awakenings -0.001 0.023 -0.046 0.045 1.002 dnorm(0,1e-3)
Morning.Feel -0.015 0.024 -0.062 0.032 1.003 dnorm(0,1e-3)
Webcm.Brghtnss -0.020 0.022 -0.062 0.023 1.008 dnorm(0,1e-3)
Webcam.Color.L -0.020 0.022 -0.062 0.023 1.009 dnorm(0,1e-3)
Webcam.Color.A -0.001 0.021 -0.044 0.04 1.000 dnorm(0,1e-3)
Webcam.Color.B -0.007 0.021 -0.049 0.034 1.001 dnorm(0,1e-3)
CloudCover 0.000 0.018 -0.037 0.036 1.000 dnorm(0,1e-3)
DayLength -0.000 0.018 -0.035 0.037 1.000 dnorm(0,1e-3)
Sleep 0.000
Light 0.000
Variances:
Post.Mean Post.SD HPD.025 HPD.975 PSRF Prior
ZQ 0.077 0.003 0.072 0.083 1.000 dgamma(1,.5)
Total.Z 0.009 0.001 0.008 0.011 1.001 dgamma(1,.5)
Time.to.Z 1.000 0.033 0.936 1.064 1.000 dgamma(1,.5)
Time.in.Wake 0.992 0.033 0.929 1.057 1.000 dgamma(1,.5)
Time.in.REM 0.473 0.016 0.444 0.505 1.000 dgamma(1,.5)
Time.in.Light 0.256 0.009 0.239 0.274 1.000 dgamma(1,.5)
Time.in.Deep 0.769 0.025 0.72 0.819 1.000 dgamma(1,.5)
Awakenings 0.933 0.031 0.872 0.992 1.000 dgamma(1,.5)
Morning.Feel 0.855 0.028 0.8 0.91 1.000 dgamma(1,.5)
Webcm.Brghtnss 0.069 0.005 0.06 0.079 1.001 dgamma(1,.5)
Webcam.Color.L 0.027 0.004 0.018 0.035 1.002 dgamma(1,.5)
Webcam.Color.A 0.999 0.030 0.941 1.059 1.000 dgamma(1,.5)
Webcam.Color.B 0.870 0.026 0.819 0.922 1.000 dgamma(1,.5)
CloudCover 0.989 0.026 0.94 1.04 1.000 dgamma(1,.5)
DayLength 0.926 0.025 0.879 0.977 1.000 dgamma(1,.5)
Sleep 0.924 0.033 0.861 0.988 1.000 dwish(iden,3)
Light 0.933 0.031 0.874 0.993 1.000 dwish(iden,3)
It takes several hours to fit using blavaan, which was getting in the way of trying to identify minimal cases, so I cut down on sampling drastically ('burnin=1000, sample=1000'). Hopefully that doesn't interfere with triggering the bug, but it made it a lot faster to iterate.
Tweaking the model repeatedly I found:
- Political/Democracy builtin dataset: 3 latent variables with 11 measurement variables: worked
- Sleep=2,Light=3,Weather=3,Music=3,IRC=2, 13: worked
- Sleep=3,Light=3,Weather=3,Music=3,IRC=2, 14: error
- Sleep+Light+Weather, with Weather reduced to 2 measurements instead of 7: error
- Sleep=9+Light=6+Weather=2, 17: error
- Sleep=8+Light=6+Weather=2, 16: error
- Sleep=7+Light=6+Weather=2, 15: worked
- Sleep=7+Light=6+Weather=3, 16: error
- Sleep=6+Light=6+Weather=3, 15: worked
- Sleep=6+Light=5+Weather=4, 15: worked
- Sleep=6+Light=5+Weather=5, 16: error
- Sleep=5+Light=5+Weather=5, 15: worked
- Sleep=4+Light=4+Weather=4+Music=3, 15: error
- Sleep=9,Weather=7, 16: worked
- Sleep=9,Weather=9, 18: error
- Sleep=8,Weather=9, 17: error
- Sleep=8,Weather=8, 16: worked
So with 5 latent variables, blavaan will fit models with <=13 measurement variables. With 3, it'll fit <=15 With 2 latent variables, it'll fit <=16 measurement variables. (I don't have any latent variables with >=17 measurement variables, so I couldn't check the singular case.) So 5+13, 3+15, or 2+16. This suggests a pattern of # of latent variables + # of measurement variables <= 18.
I looked at coeffun's call in 'bsem', but aside from indicating that it's called after JAGS is done running (successfully?), I didn't see anything that helped me:
...
verbose <- lavoptions$verbose
start.time <- proc.time()[3]
x <- NULL
if (lavm...@nx.free > 0L) {
if (!trans.exists) {
jagtrans <- try(lav2jags(model = lavpartable, lavdata = lavdata,
ov.cp = ov.cp, lv.cp = lv.cp, lv.x.wish = lavoptions$auto.cov.lv.x,
dp = dp, n.chains = n.chains, jagextra = jagextra,
inits = inits), silent = TRUE)
}
if (!inherits(jagtrans, "try-error")) {
if (jagfile) {
dir.create(path = jagdir, showWarnings = FALSE)
cat(jagtrans$model, file = paste(jagdir, "/sem.jag",
sep = ""))
save(jagtrans, file = paste(jagdir, "/semjags.rda",
sep = ""))
}
lavpartable <- mergejag(lavpartable, jagtrans$coefvec)
sampparms <- jagtrans$coefvec[, 1]
if ("monitor" %in% names(jagextra)) {
sampparms <- c(sampparms, jagextra$monitor)
}
if (inits == "jags")
jagtrans$inits <- NA
rjarg <- with(jagtrans, list(model = paste(model),
monitor = sampparms, data = data, inits = inits))
rjarg <- c(rjarg, mfj, jagcontrol)
if (suppressMessages(requireNamespace("modeest")))
runjags.options(mode.continuous = TRUE)
if (jag.do.fit) {
res <- try(do.call("run.jags", rjarg))
}
else {
res <- NULL
}
if (inherits(res, "try-error")) {
if (!trans.exists) {
dir.create(path = jagdir, showWarnings = FALSE)
cat(jagtrans$model, file = paste(jagdir, "/sem.jag",
sep = ""))
save(jagtrans, file = paste(jagdir, "/semjags.rda",
sep = ""))
}
stop("blavaan ERROR: problem with jags estimation. The jags model and data have been exported.")
}
}
else {
stop("blavaan ERROR: problem with translation from lavaan to jags.")
}
parests <- coeffun(lavpartable, res)
x <- parests$x
lavpartable <- parests$lavpartable
...
ecme...@gmail.com
Apr 25, 2016, 11:15:52 PM4/25/16
to blavaan
Gwern,
Thanks for the detailed message.
I tracked this error down to a default in the runjags package (which is used to call JAGS from R): after you reach a certain number of modelled parameters, runjags does not automatically create parameter summaries (which was causing blavaan to fail). I will have this fixed in the devel version tomorrow.
Ed
Thanks for the detailed message.
I tracked this error down to a default in the runjags package (which is used to call JAGS from R): after you reach a certain number of modelled parameters, runjags does not automatically create parameter summaries (which was causing blavaan to fail). I will have this fixed in the devel version tomorrow.
Ed
ecme...@gmail.com
Apr 26, 2016, 9:32:00 AM4/26/16
to blavaan
I believe this problem is now fixed in the devel version:
install.packages("blavaan", repos="http://faculty.missouri.edu/~merklee", type="source")
Ed
install.packages("blavaan", repos="http://faculty.missouri.edu/~merklee", type="source")
Ed
Reply all
Reply to author
Forward
0 new messages