Variable Effect _p_-value Better
----------- --------- --------- -----------------
ZQ 4.84 0.24 Yes
Total Z 19.1 0.41 Yes
Time to Z -0.65 0.69 Yes
Time in Wake -3.47 0.37 Yes
Time in REM 12.38 0.12 Yes
Time in Light 3.55 0.82 ?
Time in Deep 3.09 0.39 Yes
Awakenings -0.69 0.26 Yes

Nothing approaches the usual arbitrary cutoff, but it's worth noting that (almost) all variables are in the predicted & desired direction.
For my [CO2 sleep analysis](/zeo/CO2), I looked into more sophisticated modeling techniques and how to deal with the messiness of the Zeo data; the Zeo sleep variables are not independent of each other, have some skew, and definitely have a lot of measurement error in them, so just tossing them into a linear model isn't optimal.
After some experimentation, I defined transforms for each variable to make them normal (simply standardizing them/using `scale()` is insufficient), and tried to put together a structural equation model which would reflect the noise in the Zeo measurements (substantial and increasing over time) and also try to combine them in some hierarchical way reflecting better/worse overall sleep quality as a latent factor model.
Putting it all together:
~~~{.R}
## transform:
zma$Total.Z.2 <- zma$Total.Z^2
zma$ZQ.2 <- zma$ZQ^2
zma$Time.in.REM.2 <- zma$Time.in.REM^2
zma$Time.in.Light.2 <- zma$Time.in.Light^2
zma$Time.in.Wake.log <- log1p(zma$Time.in.Wake)
zma$Time.to.Z.log <- log1p(zma$Time.to.Z)
model1 <- '
## single-indicator measurement error model for each sleep variable assuming decent reliability:
ZQ.2_latent =~ 1*ZQ.2
ZQ.2 ~~ 0.7*ZQ.2
Total.Z.2_latent =~ 1*Total.Z.2
Total.Z.2 ~~ 0.7*Total.Z.2
Time.in.REM.2_latent =~ 1*Time.in.REM.2
Time.in.REM.2 ~~ 0.7*Time.in.REM.2
Time.in.Light.2_latent =~ 1*Time.in.Light.2
Time.in.Light.2 ~~ 0.7*Time.in.Light.2
Time.in.Deep_latent =~ 1*Time.in.Deep
Time.in.Deep ~~ 0.7*Time.in.Deep
Time.to.Z.log_latent =~ 1*Time.to.Z.log
Time.to.Z.log ~~ 0.7*Time.to.Z.log
Time.in.Wake.log_latent =~ 1*Time.in.Wake.log
Time.in.Wake.log ~~ 0.7*Time.in.Wake.log
Awakenings_latent =~ 1*Awakenings
Awakenings ~~ 0.7*Awakenings
GOOD_SLEEP =~ ZQ.2_latent + Total.Z.2_latent + Time.in.REM.2_latent + Time.in.Light.2_latent + Time.in.Deep_latent + Morning.Feel
BAD_SLEEP =~ Time.to.Z.log_latent + Time.in.Wake.log_latent + Awakenings_latent
GOOD_SLEEP ~ ZMA
BAD_SLEEP ~ ZMA
'
## Fit frequentist SEM for comparison:
library(lavaan)
l1 <- sem(model1, missing="FIML", data=scale(zma[-1])); summary(l1)
# ...Latent Variables:
# Estimate Std.Err z-value P(>|z|)
# ZQ.2_latent =~
# ZQ.2 1.000
# Total.Z.2_latent =~
# Total.Z.2 1.000
# Time.in.REM.2_latent =~
# Time.in.REM.2 1.000
# Time.in.Light.2_latent =~
# Time.in.Lght.2 1.000
# Time.in.Deep_latent =~
# Time.in.Deep 1.000
# Time.to.Z.log_latent =~
# Time.to.Z.log 1.000
# Time.in.Wake.log_latent =~
# Time.in.Wak.lg 1.000
# Awakenings_latent =~
# Awakenings 1.000
# GOOD_SLEEP =~
# ZQ.2_latent 1.000
# Total.Z.2_ltnt 1.168 0.035 33.831 0.000
# Tm.n.REM.2_ltn 0.860 0.060 14.356 0.000
# Tm.n.Lght.2_lt 0.979 0.053 18.451 0.000
# Time.n.Dp_ltnt 0.154 0.073 2.099 0.036
# Morning.Feel -0.027 0.034 -0.794 0.427
# BAD_SLEEP =~
# Tm.t.Z.lg_ltnt 1.000
# Tm.n.Wk.lg_ltn 2.321 0.520 4.464 0.000
# Awakenngs_ltnt 2.844 0.623 4.567 0.000
#
# Regressions:
# Estimate Std.Err z-value P(>|z|)
# GOOD_SLEEP ~
# ZMA -0.012 0.045 -0.264 0.792
# BAD_SLEEP ~
# ZMA -0.038 0.032 -1.199 0.231
#
# Covariances:
# Estimate Std.Err z-value P(>|z|)
# .GOOD_SLEEP ~~
# .BAD_SLEEP 0.158 0.042 3.713 0.000
## Fit Bayesian SEM:
library(blavaan)
s1 <- bsem(model1, n.chains=8, burnin=10000, sample=20000, test="none",
dp = dpriors(nu = "dnorm(0,1)", alpha = "dnorm(0,1)", beta = "dnorm(0,200)"),
jagcontrol=list(method="rjparallel"), fixed.x=FALSE, data=scale(zma[-1])); summary(s1)
# ...Latent Variables:
# Estimate Post.SD HPD.025 HPD.975 PSRF Prior
# ZQ.2_latent =~
# ZQ.2 1.000
# Total.Z.2_latent =~
# Total.Z.2 1.000
# Time.in.REM.2_latent =~
# Time.in.REM.2 1.000
# Time.in.Light.2_latent =~
# Time.in.Lght.2 1.000
# Time.in.Deep_latent =~
# Time.in.Deep 1.000
# Time.to.Z.log_latent =~
# Time.to.Z.log 1.000
# Time.in.Wake.log_latent =~
# Time.in.Wak.lg 1.000
# Awakenings_latent =~
# Awakenings 1.000
# GOOD_SLEEP =~
# ZQ.2_latent 1.000
# Total.Z.2_ltnt 1.060 0.123 0.829 1.305 1.001 dnorm(0,1e-2)
# Tm.n.REM.2_ltn 0.828 0.116 0.604 1.059 1.001 dnorm(0,1e-2)
# Tm.n.Lght.2_lt 0.953 0.120 0.726 1.195 1.001 dnorm(0,1e-2)
# Time.n.Dp_ltnt 0.769 0.112 0.56 0.996 1.001 dnorm(0,1e-2)
# Morning.Feel 0.253 0.107 0.046 0.466 1.000 dnorm(0,1e-2)
# BAD_SLEEP =~
# Tm.t.Z.lg_ltnt 1.000
# Tm.n.Wk.lg_ltn 1.550 0.342 0.935 2.24 1.003 dnorm(0,1e-2)
# Awakenngs_ltnt 1.702 0.352 1.066 2.416 1.003 dnorm(0,1e-2)
#
# Regressions:
# Estimate Post.SD HPD.025 HPD.975 PSRF Prior
# GOOD_SLEEP ~
# ZMA 0.042 0.055 -0.068 0.146 1.000 dnorm(0,200)
# BAD_SLEEP ~
# ZMA -0.033 0.041 -0.114 0.046 1.000 dnorm(0,200)
#
# Covariances:
# Estimate Post.SD HPD.025 HPD.975 PSRF Prior
# .GOOD_SLEEP ~~
# .BAD_SLEEP 0.131 0.051 0.036 0.235 1.002 dbeta(1,1)
## Plot the inferred sleep quality scores:
good <- predict(s1)[,9]
bad <- predict(s1)[,10]
qplot(zma$Date, good-bad, color=as.logical(zma$ZMA)) + stat_smooth() + theme(legend.position = "none", axis.title.x=element_blank())
~~~
![Latent sleep quality factors extracted from 9 Zeo sleep experiments and compared by ZMA supplementation status](/images/zeo/zma-sleepfactors.png)
The latent variables look sensible, with the variables loading in the expected directions.
Combining the information across sleep variables sharpens the ZMA estimates in the beneficial directions (more good sleep, less bad sleep), but the graph is unconvincing and the posterior estimates are still weak, with ~77%/80% probabilities respectively of desirable effects:
~~~{.R}
## munging the MCMC matrix/lists, appears to be variables 26,27 as of 2018-01-04, blavaan version 0.2-4:
# beta[9,11,1] -0.064641 0.042357 0.14972 0.042289 0.0548 -- 0.00045562 0.8 14466 0.035169 1.0001
posteriorGood <- s1@external$mcmcout$mcmc[1][,26][[1]]
# beta[10,11,1] -0.11429 -0.033242 0.046418 -0.033257 0.040809 -- 0.00049508 1.2 6795 0.11383 1.0003
posteriorBad <- s1@external$mcmcout$mcmc[1][,27][[1]]
mean(posteriorGood>0)
# [1] 0.77425
mean(posteriorBad<0)
# [1] 0.802
~~~
A _d_=0.04 (the factor is standardized) is a weak effect and not particularly exciting.
The higher end effect-size estimates might be worth paying for but are also not especially probable (eg _d_>=0.10 is only _P_=14%).
Ultimately, I didn't learn much about ZMA's effects from this experiment.
## Decision
Even without a formal decision analysis, given the estimated annual cost of ~\$90, ZMA is not worth taking on the basis of this self-experiment.
Should I run further experiments?
For this experiment, I was troubled by the amount of data which was NA and issues with my Zeo being reliable, possibly related to older issues involving inadequate wall socket voltage (bad wiring).
It would take more sampling than usual to get a meaningful amount of data^[Quickly testing by the crude approach of simply concatenating the data to itself, even doubling the data is not enough to provide high confidence, and only boosts _d_>0.10 from 14% to 19%.], and I still wouldn't be able to trust the results too much because I don't know what the problem is or how to fix or model it.
So while ZMA might be worth experimenting on further, I think ZMA (or other) experiments would be pointless until I can fix my Zeo data issues.