Miscellaneous usability questions
140 views
Skip to first unread message
Gwern Branwen
Jun 27, 2016, 12:46:40 AM6/27/16
to blavaan
I've continued using blavaan, switching to a different project while I
think about my productivity SEM problems, and I've noted down some
issues and questions I've run into:
1. why is blavaan so slow? It seems like even the simplest models with
3 or 4 latent variables and a few regressions can take 2 hours to fit.
Is this due to implementation issues in blavaan or is it inherent to
Bayesian estimation of latent variables?
2. is there any way to speed up blavaan aside from disabling the
posterior predictive checks which I've already done? For example, it
seems like if one could specify the initial values (as is possible in
most of the R JAGS interfaces), it would be possible to take the
post-means from earlier runs and reuse them in a later run and reduce
or skip entirely the burn-in phase, saving a fraction of the time.
The `bsem()` documentation is unclear about this: it says you can
take an equivalent `lavaan` fitted model and use its inits (but this
isn't too relevant to me since lavaan crashes on most of my datasets)
and I'm not sure what the final option listed means:
If ‘start’ is a
fitted object of class ‘lavaan’, the estimated values of the
corresponding parameters will be extracted, then perturbed in
the manner described above. If it is a model list, for
example the output of the ‘paramaterEstimates()’ [sic] function,
the values of the ‘est’ or ‘start’ or ‘ustart’ column
(whichever is found first) will be extracted.
`parameterEstimates()` is another lavaan function. But does this
imply that you can fit a model with `b1<-bsem()` and then later reuse
it with `parameterEstimates()` like
`b2<-bsem(inits=parameterEstimates(b1))`?
3. how does the `sample=10000` parameter interact with `n.chains`? Is
it 10k samples per chain (so doing 8 chains in parallel would yield
80k usable posterior samples) or is it 10k samples total (so 8 chains
would collect 1250 posterior samples each)?
This is important to me because given how
computationally-demanding blavaan is, I've been pondering renting big
Amazon EC2 instances with 64 or 128 CPU cores (only ~$4/hr), and
running that many chains in order to get a target number of samples
within a reasonable amount of wallclock time. (If it's *total*
samples, then I can simply specify `sample=100000` and `n.chains=128`
and let it run for the hour or two it needs and shut it down before
I've spent too much; if it's per-chain, then I need to set
`samples=100000/128` & `n.chains=128`.)
4. does blavaan support arithmetic expressions such as `^2` for
quadratic terms in the formula syntax like `lm()`? lavaan doesn't and
one has to create transformed variables outside of the model syntax,
which is a nuisance.
5. the `BF()` function is used to compare blavaan models to get Bayes
factors and help with model selection. The default example works, but
when I try using it on my own models, it often throws an indexing
error like:
Error in object1@test[[1]] : subscript out of bounds
Calls: BF
Typically I'm adding or removing variables, such as squared
variables for quadratic terms. The man page documentation doesn't
mention or hint at any limitations; what models can be compared using
`BF()` and how must they be evaluated by `bsem()`? Is there a required
number of posterior samples? Do the models need to have exactly the
same set of included variables? Does use of `test="none"` (to make it
run in a reasonable amount of time) or `fixed.x=FALSE` break `BF()`?
6. does blavaan support interaction terms? I am uncertain because
while throwing in 'X + Y + X*Y' into the model doesn't break it, the
output is not what I would expect.
I've been creating a SEM mediation model of a self-experiment in
which I investigate the effects of a randomized open door & bedroom
fan on a sleep quality factor via CO2 levels measured by a bedroom air
monitor (and also heat & humidity levels measured by two air
monitors); the fan and door can be expected to have main effects where
they reduce CO2 levels, but also a synergistic effect where the open
door + fan reduces CO2 much more - hence, an interaction. My current
model goes like this:
'SLEEP =~ ZQ.2 + Total.Z.2 + Time.to.Z.log +
Time.in.Wake.log + Time.in.REM.2 + Time.in.Light.2 + Time.in.Deep +
Awakenings + Morning.Feel
HUMIDITY =~ Humidity.N + Humidity.K
TEMPERATURE =~ Temperature.F.N + Temperature.F.K
CO2.N.root ~ Door.r + Fan.r + Door.r*Fan.r
HUMIDITY ~ Door.r + Fan.r + Door.r*Fan.r
TEMPERATURE ~ Door.r + Fan.r + Door.r*Fan.r
SLEEP ~ HUMIDITY + TEMPERATURE + Noise.N + CO2.N.root + Door.r +
Fan.r + Door.r*Fan.r
MP ~ CO2.N.root + SLEEP
Mnemosyne.score ~ CO2.N.root + SLEEP'
The output is not great since apparently 10k samples are not
enough with the interim data, but it looks wrong:
...
Regressions:
Post.Mean Post.SD HPD.025 HPD.975 PSRF
Prior
CO2.N.root ~
Door.r 0.048 0.026 -0.004 0.1 1.015
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
dnorm(0,1e-2)
HUMIDITY ~
Door.r 0.029 0.028 -0.028 0.084 1.012
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
TEMPERATURE ~
Door.r -0.026 0.031 -0.088 0.035 1.058
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
SLEEP ~
HUMIDIT 0.311 0.460 -0.537 1.216 1.100
dnorm(0,1e-2)
TEMPERA -0.060 0.533 -1.348 0.767 1.809
dnorm(0,1e-2)
Noise.N -0.208 0.214 -0.58 0.196 1.109
dnorm(0,1e-2)
CO2.N.r 0.036 0.225 -0.381 0.471 1.024
dnorm(0,1e-2)
Door.r -0.067 0.036 -0.138 0.004 1.130
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
...
There's `Door.r`, which is the main effect for the door variable,
and there is `Fan.r (Dr.r)` which is... the interaction term? And
there is no `Fan.r` at all. This is true of all 4 regressions using
the main+interaction terms. What's going on there?
--
gwern
http://www.gwern.net
think about my productivity SEM problems, and I've noted down some
issues and questions I've run into:
1. why is blavaan so slow? It seems like even the simplest models with
3 or 4 latent variables and a few regressions can take 2 hours to fit.
Is this due to implementation issues in blavaan or is it inherent to
Bayesian estimation of latent variables?
2. is there any way to speed up blavaan aside from disabling the
posterior predictive checks which I've already done? For example, it
seems like if one could specify the initial values (as is possible in
most of the R JAGS interfaces), it would be possible to take the
post-means from earlier runs and reuse them in a later run and reduce
or skip entirely the burn-in phase, saving a fraction of the time.
The `bsem()` documentation is unclear about this: it says you can
take an equivalent `lavaan` fitted model and use its inits (but this
isn't too relevant to me since lavaan crashes on most of my datasets)
and I'm not sure what the final option listed means:
If ‘start’ is a
fitted object of class ‘lavaan’, the estimated values of the
corresponding parameters will be extracted, then perturbed in
the manner described above. If it is a model list, for
example the output of the ‘paramaterEstimates()’ [sic] function,
the values of the ‘est’ or ‘start’ or ‘ustart’ column
(whichever is found first) will be extracted.
`parameterEstimates()` is another lavaan function. But does this
imply that you can fit a model with `b1<-bsem()` and then later reuse
it with `parameterEstimates()` like
`b2<-bsem(inits=parameterEstimates(b1))`?
3. how does the `sample=10000` parameter interact with `n.chains`? Is
it 10k samples per chain (so doing 8 chains in parallel would yield
80k usable posterior samples) or is it 10k samples total (so 8 chains
would collect 1250 posterior samples each)?
This is important to me because given how
computationally-demanding blavaan is, I've been pondering renting big
Amazon EC2 instances with 64 or 128 CPU cores (only ~$4/hr), and
running that many chains in order to get a target number of samples
within a reasonable amount of wallclock time. (If it's *total*
samples, then I can simply specify `sample=100000` and `n.chains=128`
and let it run for the hour or two it needs and shut it down before
I've spent too much; if it's per-chain, then I need to set
`samples=100000/128` & `n.chains=128`.)
4. does blavaan support arithmetic expressions such as `^2` for
quadratic terms in the formula syntax like `lm()`? lavaan doesn't and
one has to create transformed variables outside of the model syntax,
which is a nuisance.
5. the `BF()` function is used to compare blavaan models to get Bayes
factors and help with model selection. The default example works, but
when I try using it on my own models, it often throws an indexing
error like:
Error in object1@test[[1]] : subscript out of bounds
Calls: BF
Typically I'm adding or removing variables, such as squared
variables for quadratic terms. The man page documentation doesn't
mention or hint at any limitations; what models can be compared using
`BF()` and how must they be evaluated by `bsem()`? Is there a required
number of posterior samples? Do the models need to have exactly the
same set of included variables? Does use of `test="none"` (to make it
run in a reasonable amount of time) or `fixed.x=FALSE` break `BF()`?
6. does blavaan support interaction terms? I am uncertain because
while throwing in 'X + Y + X*Y' into the model doesn't break it, the
output is not what I would expect.
I've been creating a SEM mediation model of a self-experiment in
which I investigate the effects of a randomized open door & bedroom
fan on a sleep quality factor via CO2 levels measured by a bedroom air
monitor (and also heat & humidity levels measured by two air
monitors); the fan and door can be expected to have main effects where
they reduce CO2 levels, but also a synergistic effect where the open
door + fan reduces CO2 much more - hence, an interaction. My current
model goes like this:
'SLEEP =~ ZQ.2 + Total.Z.2 + Time.to.Z.log +
Time.in.Wake.log + Time.in.REM.2 + Time.in.Light.2 + Time.in.Deep +
Awakenings + Morning.Feel
HUMIDITY =~ Humidity.N + Humidity.K
TEMPERATURE =~ Temperature.F.N + Temperature.F.K
CO2.N.root ~ Door.r + Fan.r + Door.r*Fan.r
HUMIDITY ~ Door.r + Fan.r + Door.r*Fan.r
TEMPERATURE ~ Door.r + Fan.r + Door.r*Fan.r
SLEEP ~ HUMIDITY + TEMPERATURE + Noise.N + CO2.N.root + Door.r +
Fan.r + Door.r*Fan.r
MP ~ CO2.N.root + SLEEP
Mnemosyne.score ~ CO2.N.root + SLEEP'
The output is not great since apparently 10k samples are not
enough with the interim data, but it looks wrong:
...
Regressions:
Post.Mean Post.SD HPD.025 HPD.975 PSRF
Prior
CO2.N.root ~
Door.r 0.048 0.026 -0.004 0.1 1.015
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
dnorm(0,1e-2)
HUMIDITY ~
Door.r 0.029 0.028 -0.028 0.084 1.012
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
TEMPERATURE ~
Door.r -0.026 0.031 -0.088 0.035 1.058
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
SLEEP ~
HUMIDIT 0.311 0.460 -0.537 1.216 1.100
dnorm(0,1e-2)
TEMPERA -0.060 0.533 -1.348 0.767 1.809
dnorm(0,1e-2)
Noise.N -0.208 0.214 -0.58 0.196 1.109
dnorm(0,1e-2)
CO2.N.r 0.036 0.225 -0.381 0.471 1.024
dnorm(0,1e-2)
Door.r -0.067 0.036 -0.138 0.004 1.130
dnorm(0,1e-2)
Fan.r (Dr.r) -0.015 0.012 -0.038 0.008 1.024
...
There's `Door.r`, which is the main effect for the door variable,
and there is `Fan.r (Dr.r)` which is... the interaction term? And
there is no `Fan.r` at all. This is true of all 4 regressions using
the main+interaction terms. What's going on there?
--
gwern
http://www.gwern.net
Ed Merkle
Jun 27, 2016, 11:01:10 PM6/27/16
to blavaan, gw...@gwern.net
Gwern,
Here are some initial thoughts on your questions (and I will revisit again soon):
1. General-purpose samplers like JAGS or BUGS will always be somewhat slow, sacrificing speed for flexibility. However, run times on the order of hours do seem extreme. I wonder whether you are running into memory issues associated with having lots of missing data (the missing data are generally sampled as extra parameters, and this is handled behind the scenes). I'd be interested if your setup is slow on the examples from the manual, using "jagcontrol=list(method='rjparallel')" if you have multiple cores. For example, I can complete the bsem() example in about 90 seconds.
2. The "sample" argument is per chain. So sample=100 with 8 chains gives you 800 total samples.
3. In general, I think you will not have much luck with arithmetic in a blavaan model call. Behavior should be the same as lavaan here.
4. Regarding the BF() function, you may get the error if you set test="none". I should make this error more informative.
Ed
Here are some initial thoughts on your questions (and I will revisit again soon):
1. General-purpose samplers like JAGS or BUGS will always be somewhat slow, sacrificing speed for flexibility. However, run times on the order of hours do seem extreme. I wonder whether you are running into memory issues associated with having lots of missing data (the missing data are generally sampled as extra parameters, and this is handled behind the scenes). I'd be interested if your setup is slow on the examples from the manual, using "jagcontrol=list(method='rjparallel')" if you have multiple cores. For example, I can complete the bsem() example in about 90 seconds.
2. The "sample" argument is per chain. So sample=100 with 8 chains gives you 800 total samples.
3. In general, I think you will not have much luck with arithmetic in a blavaan model call. Behavior should be the same as lavaan here.
4. Regarding the BF() function, you may get the error if you set test="none". I should make this error more informative.
Ed
Ed Merkle
Jul 2, 2016, 12:17:05 PM7/2/16
to blavaan, gw...@gwern.net
Here is one additional thought about speed-up: if missing data are causing the problem, it should be possible to estimate the model using only observed data (i.e., instead of sampling the missing values as parameters, we exclude them). This would allow us to use all observed values (partial cases are allowed), so it seems similar to the "full information maximum likelihood" approach. I should be able to add this approach soon. Also, I plan to revisit the inits options.
Ed
On Sunday, June 26, 2016 at 11:46:40 PM UTC-5, Gwern Branwen wrote:
Gwern Branwen
Jul 2, 2016, 1:37:19 PM7/2/16
to Ed Merkle, blavaan
2. I checked the inits behavior and it seems my guess was wrong:
blavaan can't get inits from a fitted blavaan model because while it
doesn't crash when I try, the PSRFs don't improve noticeably and it
spits out a lot of warnings, 50+ warnings of the form:
...
3: In if (inits != "jags") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags
4: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
5: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
...
50: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
...
6. The interaction stuff hasn't changed and still leaves me perplexed
whether blavaan is genuinely estimating main+interaction effects or
has some sort of parsing issue.
On Mon, Jun 27, 2016 at 11:01 PM, Ed Merkle <ecme...@gmail.com> wrote:
> 1. General-purpose samplers like JAGS or BUGS will always be somewhat slow,
> sacrificing speed for flexibility. However, run times on the order of hours
> do seem extreme. I wonder whether you are running into memory issues
> associated with having lots of missing data (the missing data are generally
> sampled as extra parameters, and this is handled behind the scenes).
I have been benchmarking this over the past week; the missingness
varies heavily by time while the MCMC should scale with the data size,
so I've been comparing timings by adding in years/_n_ to look for
non-linear changes in runtime. There looks to be *some* but not a huge
amount, and there's not more than a 50% slowdown.
Since I have Mnemosyne data back to ~2009 (2009-05-31), sleep back to
2010 (2010-12-26) MP back to ~2012 (2012-02-16), but the Kongin only
goes back to 2014 and the Netatmo only back to May 2016, there is
enormous missingness and some variables are estimated much more
precisely than others, and may not be worth causing the slowdown in
the MCMC.
So let's drop everything after the current day and before 2016 and see
how much that speeds up 100k samples (it may even make it possible to
turn back on posterior predictive checks, which would allow model
comparison through `BF()`).
Then we can gradually extend data coverage back to see how computation
time scales with years.
Just this year:
nightShort <- night[night$Date>=as.Date("2016-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort); nrow(nightShort) / nrow(night) *
15.6 ## ~7% of the data by _n_/rows, so should take ~1.04h if
computation time scales linearly with the rows
# [1] 2678
# [1] 179
# [1] 1.060194175
system.time(s5 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort[-1]))); summary(s5)
1841 seconds or ~0.51 hours, half than expected from just _n_=179.
Adding another year of data, 2015:
nightShort2015 <- night[night$Date>=as.Date("2015-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2015); nrow(nightShort2015) /
nrow(night); (nrow(nightShort2015) / nrow(night)) * 15.6
# [1] 2678
# [1] 545
# [1] 0.2035100822
# [1] 3.174757282 ## so expected to run over 3.17h if runtime is
merely proportional to rows, but will probably run less since
missingness still isn't high back in 2015...
system.time(s6 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2015[-1]))); summary(s6)
10305 seconds or 2.86h rather than 3.17h - still better but starting
to converge on predicted.
nightShort2014 <- night[night$Date>=as.Date("2014-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2014); nrow(nightShort2014) /
nrow(night); (nrow(nightShort2014) / nrow(night)) * 15.6
# [1] 2678
# [1] 911
# [1] 0.3401792382
# [1] 5.306796117
system.time(s7 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2014[-1]))); summary(s7)
24489 seconds or 6.8h, more than 5.3h - now more than predicted
nightShort2013 <- night[night$Date>=as.Date("2013-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2013); nrow(nightShort2013) /
nrow(night); (nrow(nightShort2013) / nrow(night)) * 15.6
# [1] 2678
# [1] 1279
# [1] 0.4775952203
# [1] 7.450485437
system.time(s8 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2013[-1]))); summary(s8)
26218 seconds or 7.28h - less than predicted.
nightShort2012 <- night[night$Date>=as.Date("2012-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2012); nrow(nightShort2012) /
nrow(night); (nrow(nightShort2012) / nrow(night)) * 15.6
# [1] 2678
# [1] 1645
# [1] 0.6142643764
# [1] 9.582524272
system.time(s9 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2012[-1]))); summary(s9)
32321s or 8.97h - less than predicted.
nightShort2011 <- night[night$Date>=as.Date("2011-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2011); nrow(nightShort2011) /
nrow(night); (nrow(nightShort2011) / nrow(night)) * 15.6
# [1] 2678
# [1] 1989
# [1] 0.7427184466
# [1] 11.58640777
system.time(s10 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2011[-1])))
41317s or 11.47h, around predicted.
Presumably going back to 2010 and 2009 would yield similar results
(convergence) but that would take another two days so I'm going to end
benchmarking there.
> 2. The "sample" argument is per chain. So sample=100 with 8 chains gives you
> 800 total samples.
OK. But that doesn't explain the major slowdown I observe with
parallelism on the example model; see my other email.
> 3. In general, I think you will not have much luck with arithmetic in a
> blavaan model call. Behavior should be the same as lavaan here.
It would be nice if the blavaan formula syntax could be made as
flexible as the standard R formula syntax which allows logs, roots,
powers, and the other common transformations. I suppose it's not too
important though.
> 4. Regarding the BF() function, you may get the error if you set
> test="none". I should make this error more informative.
Yes, that seems to be it. I wasn't expecting model comparison to rely
on posterior predictive checks, though.
Actually, how does that work? I've seen PPCs expounded by Gelman as a
way to reject models but it's not one of the Bayesian model comparison
techniques I've read about before like the product-space method; I'd
always seen PPCs described as *not* giving a Bayesian comparison from
which a BF could be extracted, hence the name 'Bayesian p-value'.
--
gwern
http://www.gwern.net
blavaan can't get inits from a fitted blavaan model because while it
doesn't crash when I try, the PSRFs don't improve noticeably and it
spits out a lot of warnings, 50+ warnings of the form:
...
3: In if (inits != "jags") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags
4: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
5: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
...
50: In if (inits == "prior") { ... :
the condition has length > 1 and only the first element will be used
Calls: system.time ... tryCatchList -> tryCatchOne ->
doTryCatch -> lav2jags -> set_inits
...
6. The interaction stuff hasn't changed and still leaves me perplexed
whether blavaan is genuinely estimating main+interaction effects or
has some sort of parsing issue.
On Mon, Jun 27, 2016 at 11:01 PM, Ed Merkle <ecme...@gmail.com> wrote:
> 1. General-purpose samplers like JAGS or BUGS will always be somewhat slow,
> sacrificing speed for flexibility. However, run times on the order of hours
> do seem extreme. I wonder whether you are running into memory issues
> associated with having lots of missing data (the missing data are generally
> sampled as extra parameters, and this is handled behind the scenes).
varies heavily by time while the MCMC should scale with the data size,
so I've been comparing timings by adding in years/_n_ to look for
non-linear changes in runtime. There looks to be *some* but not a huge
amount, and there's not more than a 50% slowdown.
Since I have Mnemosyne data back to ~2009 (2009-05-31), sleep back to
2010 (2010-12-26) MP back to ~2012 (2012-02-16), but the Kongin only
goes back to 2014 and the Netatmo only back to May 2016, there is
enormous missingness and some variables are estimated much more
precisely than others, and may not be worth causing the slowdown in
the MCMC.
So let's drop everything after the current day and before 2016 and see
how much that speeds up 100k samples (it may even make it possible to
turn back on posterior predictive checks, which would allow model
comparison through `BF()`).
Then we can gradually extend data coverage back to see how computation
time scales with years.
Just this year:
nightShort <- night[night$Date>=as.Date("2016-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort); nrow(nightShort) / nrow(night) *
15.6 ## ~7% of the data by _n_/rows, so should take ~1.04h if
computation time scales linearly with the rows
# [1] 2678
# [1] 179
# [1] 1.060194175
system.time(s5 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort[-1]))); summary(s5)
1841 seconds or ~0.51 hours, half than expected from just _n_=179.
Adding another year of data, 2015:
nightShort2015 <- night[night$Date>=as.Date("2015-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2015); nrow(nightShort2015) /
nrow(night); (nrow(nightShort2015) / nrow(night)) * 15.6
# [1] 2678
# [1] 545
# [1] 0.2035100822
# [1] 3.174757282 ## so expected to run over 3.17h if runtime is
merely proportional to rows, but will probably run less since
missingness still isn't high back in 2015...
system.time(s6 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2015[-1]))); summary(s6)
10305 seconds or 2.86h rather than 3.17h - still better but starting
to converge on predicted.
nightShort2014 <- night[night$Date>=as.Date("2014-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2014); nrow(nightShort2014) /
nrow(night); (nrow(nightShort2014) / nrow(night)) * 15.6
# [1] 2678
# [1] 911
# [1] 0.3401792382
# [1] 5.306796117
system.time(s7 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2014[-1]))); summary(s7)
24489 seconds or 6.8h, more than 5.3h - now more than predicted
nightShort2013 <- night[night$Date>=as.Date("2013-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2013); nrow(nightShort2013) /
nrow(night); (nrow(nightShort2013) / nrow(night)) * 15.6
# [1] 2678
# [1] 1279
# [1] 0.4775952203
# [1] 7.450485437
system.time(s8 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2013[-1]))); summary(s8)
26218 seconds or 7.28h - less than predicted.
nightShort2012 <- night[night$Date>=as.Date("2012-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2012); nrow(nightShort2012) /
nrow(night); (nrow(nightShort2012) / nrow(night)) * 15.6
# [1] 2678
# [1] 1645
# [1] 0.6142643764
# [1] 9.582524272
system.time(s9 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2012[-1]))); summary(s9)
32321s or 8.97h - less than predicted.
nightShort2011 <- night[night$Date>=as.Date("2011-01-01") &
night$Date<Sys.Date(),]
nrow(night); nrow(nightShort2011); nrow(nightShort2011) /
nrow(night); (nrow(nightShort2011) / nrow(night)) * 15.6
# [1] 2678
# [1] 1989
# [1] 0.7427184466
# [1] 11.58640777
system.time(s10 <- bsem(model3, test="none", n.chains=6,
jagcontrol=list(method="rjparallel"), sample=100000, fixed.x=FALSE,
data=scale(nightShort2011[-1])))
41317s or 11.47h, around predicted.
Presumably going back to 2010 and 2009 would yield similar results
(convergence) but that would take another two days so I'm going to end
benchmarking there.
> 2. The "sample" argument is per chain. So sample=100 with 8 chains gives you
> 800 total samples.
parallelism on the example model; see my other email.
> 3. In general, I think you will not have much luck with arithmetic in a
> blavaan model call. Behavior should be the same as lavaan here.
flexible as the standard R formula syntax which allows logs, roots,
powers, and the other common transformations. I suppose it's not too
important though.
> 4. Regarding the BF() function, you may get the error if you set
> test="none". I should make this error more informative.
on posterior predictive checks, though.
Actually, how does that work? I've seen PPCs expounded by Gelman as a
way to reject models but it's not one of the Bayesian model comparison
techniques I've read about before like the product-space method; I'd
always seen PPCs described as *not* giving a Bayesian comparison from
which a BF could be extracted, hence the name 'Bayesian p-value'.
--
gwern
http://www.gwern.net
Gwern Branwen
Jul 2, 2016, 1:53:47 PM7/2/16
to Ed Merkle, blavaan
On Sat, Jul 2, 2016 at 1:37 PM, Gwern Branwen <gw...@gwern.net> wrote:
> Since I have Mnemosyne data back to ~2009 (2009-05-31), sleep back to
> 2010 (2010-12-26) MP back to ~2012 (2012-02-16), but the Kongin only
> goes back to 2014 and the Netatmo only back to May 2016, there is
> enormous missingness and some variables are estimated much more
> precisely than others, and may not be worth causing the slowdown in
> the MCMC.
> So let's drop everything after the current day and before 2016 and see
> how much that speeds up 100k samples (it may even make it possible to
> turn back on posterior predictive checks, which would allow model
> comparison through `BF()`).
> Then we can gradually extend data coverage back to see how computation
> time scales with years.
If anyone wants to experiment with this dataset themselves, 'night' is
> Since I have Mnemosyne data back to ~2009 (2009-05-31), sleep back to
> 2010 (2010-12-26) MP back to ~2012 (2012-02-16), but the Kongin only
> goes back to 2014 and the Netatmo only back to May 2016, there is
> enormous missingness and some variables are estimated much more
> precisely than others, and may not be worth causing the slowdown in
> the MCMC.
> So let's drop everything after the current day and before 2016 and see
> how much that speeds up 100k samples (it may even make it possible to
> turn back on posterior predictive checks, which would allow model
> comparison through `BF()`).
> Then we can gradually extend data coverage back to see how computation
> time scales with years.
available at https://dl.dropboxusercontent.com/u/85192141/2016-07-02-night.csv
So this should work:
night <- read.csv("2016-07-02-night.csv")
library(blavaan)
model3 <- '
SLEEP =~ ZQ.2 + Total.Z.2 + Time.to.Z.log +
Time.in.Wake.log + Time.in.REM.2 + Time.in.Light.2 + Time.in.Deep +
Awakenings + Morning.Feel
HUMIDITY =~ Humidity.N + Humidity.K
TEMPERATURE =~ Temperature.F.N + Temperature.F.K
CO2.N.root ~ Door.r + Fan.r + Door.r*Fan.r
HUMIDITY ~ Door.r + Fan.r + Door.r*Fan.r
TEMPERATURE ~ Door.r + Fan.r + Door.r*Fan.r
SLEEP ~ HUMIDITY + TEMPERATURE + Noise.N + CO2.N.root + Door.r
+ Fan.r + Door.r*Fan.r
MP ~ CO2.N.root + SLEEP
Mnemosyne.score ~ CO2.N.root + SLEEP
'
system.time(s3 <- bsem(model3, test="none", n.chains=6,
Time.in.Wake.log + Time.in.REM.2 + Time.in.Light.2 + Time.in.Deep +
Awakenings + Morning.Feel
HUMIDITY =~ Humidity.N + Humidity.K
TEMPERATURE =~ Temperature.F.N + Temperature.F.K
CO2.N.root ~ Door.r + Fan.r + Door.r*Fan.r
HUMIDITY ~ Door.r + Fan.r + Door.r*Fan.r
TEMPERATURE ~ Door.r + Fan.r + Door.r*Fan.r
SLEEP ~ HUMIDITY + TEMPERATURE + Noise.N + CO2.N.root + Door.r
+ Fan.r + Door.r*Fan.r
MP ~ CO2.N.root + SLEEP
Mnemosyne.score ~ CO2.N.root + SLEEP
'
jagcontrol=list(method="rjparallel"), fixed.x=FALSE,
data=scale(night[-1])))
--
gwern
http://www.gwern.net
Reply all
Reply to author
Forward
0 new messages