CO2/ventilation sleep experiment - Gwern.net


Paul Christiano asked about some curious psychology results suggesting that ambient levels of CO2 considered safe and having no effects on cognition may actually have substantial harmful effects, particularly Satish et al 2012/Allen et al 2015. I looked into the research literature for related papers, but I had a hard time believing the effect: it’s hard to believe that CO2, which makes up <1% of normal air, could have such effects, especially since such effects would mean we should see large correlations with weather, altitude, test location etc, and the US Navy sponsored extensive testing of CO2 levels at levels reaching orders of magnitude larger than Satish/Allen (summary) and found cognitive effects only at very high PPMs (>25,000 PPM).

But after looking at consumer air quality sensors, I decided to get a Netatmo weather station device. It cost ~$140 and records air temperature, humidity, CO2, & barometric pressure.

The last observation interested me the most.

I previously kept my bedroom tightly sealed to block light and to reduce how much hot/cold air infiltrates during summer/winter. In the summer, I don’t want the (recirculating) AC on while I sleep because the moving air bothers me; in winter, cold air isn’t so bad, but it’s hard to get the electric wall heating units to keep the room consistently warm and I’ve taken to relying on electric blankets.

So it wouldn’t be surprising if CO2 built up in the bedroom while I slept and then when I opened it in the morning and the air mixed with the rest of the apartment, the average CO2 level increased. But if one room can mix with and increase the Netatmo’s room all the way up to 1934PPM, what must the original concentration have been like? It would have had to have been at least twice as high (~3000PPM) in which case the bedroom CO2 levels could easily have been >7x the normal levels and is starting to near levels which are considered unacceptable (5000PPM). I don’t know if CO2 has any effect on my sleep, but if it does, differences like 7x are enough that the effect could be noticeable.

High levels CO2 could potentially worsen sleep, similar to hypercapnia, where the panic/arousal impedes sleep. Existing sleep research mostly investigates sleep apnea, infants, and SIDS.

First, to get an idea of normal CO2 levels during sleep, I moved the Netatmo to my bedroom for one night and turned everything off & shut the room like usual. CO2 climbed gradually to 2832PPM at 5:51AM and remained elevated until I woke up and opened the door. The second night, I did the same but left the door open. CO2 peaked at 1619PPM at 5:46AM and likewise only fell when the fans/AC/dehumidifier turned on in the morning. For the third night, door open, but one fan turned on slightly facing the door. CO2 peaked at 1405PPM at 4:58AM. The fourth night reverted to closed, first peaking at 2775PPM at 5:21AM and then peaking at 2912PPM at 9:55AM. Fifth night was open with fan: 1407PPM at 5:57AM. When the door is closed and the fan is on, there is a small but not useful reduction in CO2 (peaking one night at 2452PPM), which could be due either to a more even distribution within the room or possibly from accelerating air infiltration from the outside. The next 5 nights yielded similar numbers: closed leads to steady increases peaking 2900-3000PPM, open with no fan or with fan tends to not be more than ~1400PPM and possibly less.

To model sleep effects, I want to take into account bedroom temperature, humidity, CO2, sound, and measure sleep. The Netatmo measures temperature/humidity/CO2/sound. My Kongin temperature/humidity logger, which I have used since October 2014, will be useful for examining temperature/humidity effects stretching back much further than the Netatmo’s measurements and can be used to double-check the Netatmo measurements. For sleep, I have the Zeo.

The door and fan cannot be blinded since I do not have access to a super-quiet fan which could possibly turn on at random while I’m asleep, so this would be a non-blind randomized experiment. Since the windows are sealed, and the door and fan can either be on or not, that leads to 4 possible conditions: door closed fan off, door closed fan on, door open fan off, door open fan on. So randomization can be done in blocks of 4, covering each possible combination, and recorded as binary variables. I don’t want to omit the fan-only condition since that helps check for possible harm to my sleep from the fan’s noise/air-movement.

Power/sample size is unclear since I haven’t been able to get most of the papers reporting on CO2/sleep data, but the most recent paper, Strøm-Tejsen et al 2014b, while not reporting any real effect sizes, claims a number of statistically-significant effects with 16 subjects with what seems to be less than a month of data each (and possibly as little as 8 days total); 3 months is probably adequate.

Changing CO2 levels may have direct effects on sleep, but it may also be that the change in temperature/humidity is responsible, or that having the door open or the fan running are themselves directly responsible and any effect of CO2 is merely because opening the door does multiple things. This sounds like a mediation model, in which we want to see if the intervention has an effect on sleep quality only when mediated through an intermediate variable, in this case, CO2. So writing it down in Lavaan syntax, something like:

The regression CO2 ~ Door.r + Fan.r + Door.r*Fan.r will confirm that the randomized intervention reduced CO2 levels by a certain amount (at least 1000PPM in our case); the exact reduction will vary for each night, while the intervention is exactly the same, so which is a better predictor will indicate which is more important causally and help distinguish between scenarios. (For example, it could be that the intervention of leaving the door open and running a fan has bad effects in moving air around enough to disturb me, but the lower CO2 the better; in that case, the coefficients of CO2 ~ Door.r + Fan.r + Door.r*Fan.r will be highly negative, Sleep ~ CO2 highly positive, and Sleep ~ Door.r + Fan.r + Door.r*Fan.r will be highly negative. Whether I would want to start leaving the door open would then depend on whether CO2 ~ Door.r + Fan.r + Door.r*Fan.r * Sleep ~ CO2 > Sleep ~ Door.r + Fan.r + Door.r*Fan.r - that is, the net improvement on sleep through the indirect path of reduced CO2 and CO2 to sleep is better than the worsening through the direct path of intervention to sleep. We could also skip the CO2 variable entirely and model it as a simple two-group comparison of means, Sleep ~ Door.r + Fan.r + Door.r*Fan.r, but then we wouldn’t know what role CO2 played: if CO2 really does worsen sleep but the fan destroys the benefits, then we can think about alternatives ways to reduce CO2 without using a fan, like removing window coverings so I can open windows at night.)

The inclusion of average nightly humidity & temperature may not be useful; my regressions beforehand did not turn up any linear or quadratic relationship between humidity/temperature and the sleep latent variable, ZQ, Total.Z, or Morning.Feel, and plotting did not show any obvious relationship.

Netatmo accuracy - indoor heat: Accuracy: ± 0.3°C / ± 0.54°F - indoor humidity: Accuracy: ± 3% - indoor CO2: Accuracy: ± 50 ppm or ± 5%

Netatmo usage started: 28 May 2016

zeo <- read.csv("~/wiki/docs/zeo/gwern-zeodata.csv")
zeo$Date <- as.Date(sapply(strsplit(as.character(zeo$Rise.Time), " "), function(x) { x[1] }), format="%m/%d/%Y")
zeo$Rise.Time <- sapply(strsplit(as.character(zeo$Rise.Time), " "), function(x) { x[2] })
zeo$Start.of.Night <- sapply(strsplit(as.character(zeo$Start.of.Night), " "), function(x) { x[2] })
interval <- function(x) { if (!is.na(x)) { if (grepl(" s",x)) as.integer(sub(" s","",x)) else { y <- unlist(strsplit(x, ":")); as.integer(y[[1]])*60 + as.integer(y[[2]]); } } else NA }
zeo$Rise.Time <- sapply(zeo$Rise.Time, interval)
zeo$Start.of.Night <- sapply(zeo$Start.of.Night, interval)
zeo <- zeo[!is.na(zeo$Date),]
zeo[(zeo$Date >= as.Date("2013-03-11")),]$Rise.Time <- (zeo[(zeo$Date >= as.Date("2013-03-11")),]$Rise.Time + 226) %% (24*60)
zeo[(zeo$Date >= as.Date("2013-03-11")),]$Start.of.Night <- (zeo[(zeo$Date >= as.Date("2013-03-11")),]$Start.of.Night + 226+680)
zeo$Start.of.Night <- ifelse(zeo$Start.of.Night<500, zeo$Start.of.Night + 1440, zeo$Start.of.Night) # undo wrapping around to 0
zeo$Time.to.Z.log <- log1p(zeo$Time.to.Z) # reduce skew
zeo <- subset(zeo, select=c("Date", "ZQ", "Total.Z", "Time.to.Z.log", "Time.in.Wake", "Time.in.REM", "Time.in.Light", "Time.in.Deep", "Awakenings", "Morning.Feel")) ## NOTE: need to create a fresh `humidity-all.csv` and `netatmo/all.csv` using existing R routines *before* running this:
kongin <- read.csv("~/selfexperiment/humidity/humidity-all.csv", colClasses=c("POSIXct", "numeric", "numeric"))
kongin$Date <- as.Date(kongin$Timestamp, tz="EST")
## only interested in midnight-to-10AM:
kongin <- kongin[as.integer(strftime(kongin$Timestamp, format="%H")) < 10,]
## aggregate data both by average and by maximum; peak temperatures/humidity may damage sleep as much or more than high averages
konginDailyMean <- aggregate(cbind(Humidity, Temperature.C) ~ Date, mean, data=kongin)
colnames(konginDailyMean) <- c("Date", "Humidity.K.mean", "Temperature.K.mean")
konginDailyMax <- aggregate(cbind(Humidity, Temperature.C) ~ Date, max, data=kongin)
colnames(konginDailyMax) <- c("Date", "Humidity.K.max", "Temperature.K.max")
konginDaily <- merge(konginDailyMean, konginDailyMax, all=TRUE)
summary(konginDaily) netatmo <- read.csv("~/selfexperiment/netatmo/all.csv")
netatmo$Date <- as.Date(netatmo$Date)
netatmo <- netatmo[as.integer(strftime(netatmo$Timestamp, format="%H")) < 10,]
netatmo <- netatmo[netatmo$Date >= as.Date("2016-06-02"),] # include only experiment data from when Netatmo was in bedroom at night
netatmoDailyMean <- aggregate(cbind(Humidity, Temperature, CO2, Noise) ~ Date, mean, data=netatmo)
colnames(netatmoDailyMean) <- c("Date", "Humidity.N.mean", "Temperature.N.mean", "CO2.N.mean", "Noise.N.mean")
netatmoDailyMax <- aggregate(cbind(Humidity, Temperature, CO2, Noise) ~ Date, max, data=netatmo)
colnames(netatmoDailyMax) <- c("Date", "Humidity.N.max", "Temperature.N.max", "CO2.N.max", "Noise.N.max")
netatmoDaily <- merge(netatmoDailyMean, netatmoDailyMax, all=TRUE) airNight <- merge(konginDaily, netatmoDaily, all=TRUE) door <- read.csv(stdin(), header=TRUE, colClasses=c("Date","integer","integer", "integer"))
Date,Fan.r,Door.r,Mattress.r
2016-06-03,0,0,NA
2016-06-04,0,1,NA
2016-06-05,1,1,NA
2016-06-06,0,0,NA
2016-06-07,1,1,NA
2016-06-08,0,1,NA
2016-06-09,NA,NA,NA
2016-06-10,1,1,NA
2016-06-11,NA,NA,NA
2016-06-12,0,0,NA
2016-06-13,0,1,NA
2016-06-14,NA,NA,NA
2016-06-15,1,0,NA
2016-06-16,0,0,NA
2016-06-17,1,0,NA
2016-06-18,0,1,NA
2016-06-19,1,1,NA
2016-06-20,NA,NA,NA
2016-06-21,1,0,NA
2016-06-22,0,0,NA
2016-06-23,0,1,NA
2016-06-24,1,1,NA
2016-06-25,1,1,NA
2016-06-26,0,0,NA
2016-06-27,0,0,NA
2016-06-28,1,0,NA
2016-06-29,0,1,NA
2016-06-30,1,1,NA
2016-07-01,1,0,NA
2016-07-02,0,1,NA
2016-07-03,0,0,NA
2016-07-04,0,0,NA
2016-07-05,1,1,NA
2016-07-06,1,0,NA
2016-07-07,0,1,NA
2016-07-08,1,0,NA
2016-07-09,0,1,NA
2016-07-10,1,1,NA
2016-07-11,0,0,NA
2016-07-12,1,1,NA
2016-07-13,0,0,NA
2016-07-14,1,0,NA
2016-07-15,0,1,NA
2016-07-16,0,1,NA
2016-07-17,0,0,NA
2016-07-18,1,1,NA
2016-07-19,1,0,NA
2016-07-20,NA,NA,NA
2016-07-21,1,0,NA
2016-07-22,1,1,NA
2016-07-23,0,0,NA
2016-07-24,0,1,NA
2016-07-25,1,1,NA
2016-07-26,1,0,NA
2016-07-27,0,1,NA
2016-07-28,0,0,NA
2016-07-29,NA,NA,NA
2016-07-30,0,0,NA
2016-07-31,1,0,NA
2016-08-01,1,1,NA
2016-08-02,0,1,NA
2016-08-03,0,1,NA
2016-08-04,1,0,NA
2016-08-05,0,0,NA
2016-08-06,1,1,NA
2016-08-07,0,1,NA
2016-08-08,1,1,NA
2016-08-09,1,0,NA
2016-08-10,0,0,NA
2016-08-11,1,0,NA
2016-08-12,0,1,NA
2016-08-13,0,0,NA
2016-08-14,1,1,NA
2016-08-15,0,0,NA
2016-08-16,1,1,NA
2016-08-17,1,0,NA
2016-08-18,0,1,NA
2016-08-19,0,1,NA
2016-08-20,1,0,NA
2016-08-21,1,1,NA
2016-08-22,0,0,NA
2016-08-23,0,0,NA
2016-08-24,1,0,NA
2016-08-25,0,1,NA
2016-08-26,1,1,NA
2016-08-27,0,1,NA
2016-08-28,1,0,NA
2016-08-29,1,1,NA
2016-08-30,0,0,NA
2016-08-31,0,0,NA
2016-09-01,1,1,NA
2016-09-02,0,1,NA
2016-09-03,1,0,NA
2016-09-04,0,1,NA
2016-09-05,1,1,NA
2016-09-06,1,0,NA
2016-09-07,0,0,NA
2016-09-08,0,1,NA
2016-09-09,1,1,NA
2016-09-10,1,0,NA
2016-09-11,0,0,NA
2016-09-12,0,1,NA
2016-09-13,1,1,NA
2016-09-14,0,0,NA
2016-09-15,1,0,NA
2016-09-16,1,0,NA
2016-09-17,NA,NA,NA
2016-09-18,0,0,NA
2016-09-19,1,1,NA
2016-09-20,0,0,NA
2016-09-21,0,1,NA
2016-09-22,1,1,NA
2016-09-23,0,0,NA
2016-09-24,0,1,NA
2016-09-25,1,0,NA
2016-09-26,0,1,NA
2016-09-27,0,0,NA
2016-09-28,0,1,NA
2016-09-29,1,0,NA
2016-09-30,1,1,NA
2016-10-01,0,0,NA
2016-10-02,1,1,NA
2016-10-03,0,1,NA
2016-10-04,1,1,NA
2016-10-05,0,1,NA
2016-10-06,0,1,NA
2016-10-07,1,0,NA
2016-10-08,0,0,NA
2016-10-09,1,1,NA
2016-10-10,1,0,NA
2016-10-11,0,0,NA
2016-10-12,0,1,NA
2016-10-13,1,1,NA
2016-10-14,0,0,NA
2016-10-15,1,1,NA
2016-10-16,0,1,NA
2016-10-17,1,0,NA
2016-10-18,1,1,NA
2016-10-19,1,0,NA
2016-10-20,0,1,NA
2016-10-21,1,1,NA
2016-10-22,0,0,NA
2016-10-23,0,0,NA
2016-10-24,0,1,NA
2016-10-25,1,1,NA
2016-10-26,1,0,NA
2016-10-27,0,1,NA
2016-10-28,1,1,0
2016-10-29,1,0,1
2016-10-30,0,0,1
2016-10-31,0,0,0
2016-11-01,1,0,1
2016-11-02,0,1,0
2016-11-03,1,1,0
2016-11-04,0,1,1
2016-11-05,1,1,1
2016-11-06,0,0,0
2016-11-07,1,0,1
2016-11-08,1,0,0
2016-11-09,0,1,0
2016-11-10,1,0,1
2016-11-11,NA,NA,0
2016-11-12,0,0,1
2016-11-13,1,0,0
2016-11-14,1,1,1
2016-11-15,0,1,1
2016-11-16,1,0,0
2016-11-17,1,0,1
2016-11-18,0,1,0
2016-11-19,1,1,0
2016-11-20,0,0,1
2016-11-21,0,1,0
2016-11-22,1,0,1
2016-11-23,1,1,1
2016-11-24,0,0,0
2016-11-25,0,0,0
2016-11-26,1,1,1
2016-11-27,0,1,0
2016-11-28,0,1,1
2016-11-29,0,1,0
2016-11-30,1,1,0
2016-12-01,0,0,1
2016-12-02,1,1,1
2016-12-03,0,0,0
2016-12-04,0,1,0
2016-12-05,1,1,1
2016-12-06,1,0,1
2016-12-07,0,1,0
2016-12-08,0,0,1
2016-12-09,1,1,0
2016-12-10,NA,NA,NA
2016-12-11,NA,NA,NA
2016-12-12,NA,NA,NA
2016-12-13,NA,NA,NA
2016-12-14,NA,NA,NA
2016-12-15,NA,NA,NA
2016-12-16,NA,NA,NA
2016-12-17,NA,NA,NA
2016-12-18,NA,NA,NA
2016-12-19,NA,NA,NA
2016-12-20,NA,NA,NA
2016-12-21,NA,NA,NA
2016-12-22,NA,NA,NA
2016-12-23,NA,NA,NA
2016-12-24,NA,NA,NA
2016-12-25,NA,NA,NA
2016-12-26,NA,NA,NA
2016-12-27,NA,NA,NA
2016-12-28,NA,NA,NA
2016-12-29,NA,NA,NA
2016-12-30,NA,NA,NA
2016-12-31,1,1,0
2017-01-01,NA,NA,NA
2017-01-02,0,1,1
2017-01-03,1,0,1
2017-01-04,0,0,0
2017-01-05,1,1,0
2017-01-06,1,0,0
2017-01-07,0,0,1
2017-01-08,0,0,0
2017-01-09,NA,NA,1
2017-01-10,0,1,1
2017-01-11,0,0,0
2017-01-12,0,0,0
2017-01-13,1,1,1
2017-01-14,1,0,1
2017-01-15,0,0,0
2017-01-16,0,1,0
2017-01-17,1,1,1
2017-01-18,1,1,1
2017-01-19,0,1,0
2017-01-20,0,1,0
2017-01-21,1,0,1
2017-01-22,0,1,0
2017-01-23,1,1,0
2017-01-24,0,0,1
2017-01-25,1,0,1
2017-01-26,1,0,0
2017-01-27,1,1,0
2017-01-28,0,0,1
2017-01-29,0,1,1
2017-01-30,1,1,0
2017-01-31,0,0,0
2017-02-01,0,1,1
2017-02-02,1,0,1
2017-02-03,0,1,0
2017-02-04,0,0,0
2017-02-05,1,0,0
2017-02-06,1,1,1
2017-02-07,1,1,1
2017-02-08,1,0,1
2017-02-09,1,0,1
2017-02-10,0,1,1
2017-02-11,1,1,NA
2017-02-12,0,0,NA
2017-02-13,0,0,NA
2017-02-14,1,1,NA
2017-02-15,1,0,NA
2017-02-16,0,1,NA
2017-02-17,1,1,NA
2017-02-18,0,0,NA
2017-02-19,0,1,NA
2017-02-20,1,0,NA
2017-02-21,0,1,NA
2017-02-22,1,0,NA
2017-02-23,1,1,NA
2017-02-24,0,0,NA
2017-02-25,0,0,NA
2017-02-26,1,1,NA
2017-02-27,NA,NA,NA
2017-02-28,1,0,NA
2017-03-01,0,0,NA
2017-03-02,1,1,NA
2017-03-03,0,0,NA
2017-03-04,NA,NA,NA
2017-03-05,0,1,NA
2017-03-06,1,0,NA
2017-03-07,0,1,NA
2017-03-08,1,1,NA
2017-03-09,0,0,NA
2017-03-10,0,0,NA
2017-03-11,0,1,NA
2017-03-12,1,1,NA
2017-03-13,0,0,NA
2017-03-14,1,0,NA
2017-03-15,0,1,NA
2017-03-16,0,0,NA
2017-03-17,1,1,NA
2017-03-18,1,0,NA
2017-03-19,0,0,NA
2017-03-20,0,1,NA
2017-03-21,1,0,NA
2017-03-22,1,1,NA
2017-03-23.1,0,NA
2017-03-24,0,1,NA
2017-03-25,0,0,NA
2017-03-26,1,1,NA
2017-03-27,0,1,NA
2017-03-28,0,0,NA
2017-03-29,1,0,NA
2017-03-30,1,0,NA
2017-03-31,0,1,NA
2017-04-01,0,0,NA
2017-04-02,1,1,NA
2017-04-03,1,1,NA
2017-04-04,0,1,NA
2017-04-05,0,0,NA
2017-04-06,1,0,NA
2017-04-07,1,1,NA
2017-04-08,0,1,NA
2017-04-09,1,0,NA
2017-04-10,0,0,NA
2017-04-11,1,0,NA
2017-04-12,0,0,NA
2017-04-13,NA,NA,NA
2017-04-14,NA,NA,NA
2017-04-15,NA,NA,NA
2017-04-16,NA,NA,NA
2017-04-17,NA,NA,NA
2017-04-18,NA,NA,NA
2017-04-19,0,0,NA
2017-04-20,0,1,NA
2017-04-21,1,0,NA
2017-04-22,1,1,NA
2017-04-23,0,1,NA
2017-04-24,0,0,NA
2017-04-25,1,1,NA
2017-04-26,1,0,NA
2017-04-27,1,1,NA
2017-04-28,1,1,NA
2017-04-29,1,0,NA
2017-04-30,NA,NA,NA
2017-05-01,0,1,NA
2017-05-02,0,0,NA
2017-05-03,1,0,NA
2017-05-04,0,1,NA
2017-05-05,1,1,NA
2017-05-06,0,0,NA
2017-05-07,0,0,NA
2017-05-08,1,1,NA
2017-05-09,0,1,NA
2017-05-10,1,0,NA
2017-05-11,0,0,NA
2017-05-12,0,1,NA
2017-05-13,1,0,NA
2017-05-14,1,1,NA
2017-05-15,1,1,NA
2017-05-16,0,0,NA
2017-05-17,1,0,NA
2017-05-18,0,1,NA
2017-05-19,0,1,NA
2017-05-20,1,1,NA
2017-05-21,0,1,NA
2017-05-22,0,1,NA
2017-05-23,1,1,NA
2017-05-24,0,1,NA
2017-05-25,0,0,NA
2017-05-26,1,1,NA
2017-05-27,0,1,NA
2017-05-28,1,1,NA
2017-05-29,0,1,NA
2017-05-30,0,0,NA
2017-05-31,0,1,NA
2017-06-01,1,1,NA
2017-06-02,0,1,NA
2017-06-03,NA,NA,NA
2017-06-04,NA,NA,NA
2017-06-05,NA,NA,NA
2017-06-06,NA,NA,NA
2017-06-07,NA,NA,NA
2017-06-08,NA,NA,NA
2017-06-09,NA,NA,NA
2017-06-10,NA,NA,NA
2017-06-11,NA,NA,NA
2017-06-12,NA,NA,NA
2017-06-13,NA,NA,NA
2017-06-14,NA,NA,NA
2017-06-15,NA,NA,NA
2017-06-16,NA,NA,NA
2017-06-17,NA,NA,NA
2017-06-18,NA,NA,NA
2017-06-19,NA,NA,NA
2017-06-20,NA,NA,NA
2017-06-21,NA,NA,NA
2017-06-22,NA,NA,NA
2017-06-23,NA,NA,NA
2017-06-24,NA,NA,NA
2017-06-25,NA,NA,NA
2017-06-26,NA,NA,NA
2017-06-27,NA,NA,NA
2017-06-28,NA,NA,NA
2017-06-29,NA,NA,NA
2017-06-30,NA,NA,NA
2017-07-01,NA,NA,NA
2017-07-02,NA,NA,NA
2017-07-03,NA,NA,NA
2017-07-04,NA,NA,NA
2017-07-05,NA,NA,NA
2017-07-06,NA,NA,NA
2017-07-07,NA,NA,NA
2017-07-08,NA,NA,NA
2017-07-09,NA,NA,NA
2017-07-10,NA,NA,NA
2017-07-11,NA,NA,NA
2017-07-12,0,0,NA
2017-07-13,1,0,NA
2017-07-14,1,1,NA
2017-07-15,0,1,NA
2017-07-16,0,1,NA
2017-07-17,0,0,NA
2017-07-18,1,0,NA
2017-07-19,1,1,NA
2017-07-20,0,0,NA
2017-07-21,1,0,NA
2017-07-22,0,1,NA
2017-07-23,1,1,NA
2017-07-24,0,0,NA
2017-07-25,1,0,NA
2017-07-26,0,0,NA
2017-07-27,0,1,NA
2017-07-28,1,0,NA
2017-07-29,1,1,NA
2017-07-30,1,0,NA
2017-07-31,1,1,NA
2017-08-01,0,0,NA
2017-08-02,0,1,NA
2017-08-03,1,0,NA
2017-08-04,0,0,NA
2017-08-05,1,1,NA
2017-08-06,0,1,NA
2017-08-07,1,0,NA
2017-08-08,0,0,NA
2017-08-09,0,1,NA
2017-08-10,1,1,NA
2017-08-11,0,1,NA
2017-08-12,1,1,NA
2017-08-13,0,0,NA
2017-08-14,1,0,NA
2017-08-15,0,0,NA
2017-08-16,1,0,NA
2017-08-17,1,1,NA
2017-08-18,0,1,NA
2017-08-19,0,1,NA
2017-08-20,1,1,NA
2017-08-21,0,0,NA
2017-08-22,1,0,NA
2017-08-23,1,0,NA
2017-08-24,0,1,NA
2017-08-25,1,1,NA
2017-08-26,0,0,NA
2017-08-27,0,1,NA
2017-08-28,0,0,NA
2017-08-29,0,1,NA
2017-08-30,0,1,NA
2017-08-31,1,1,NA
2017-09-01,0,1,NA
2017-09-02,1,1,NA
2017-09-03,0,0,NA
2017-09-04,1,0,NA
2017-09-05,0,0,NA
2017-09-06,1,1,NA
2017-09-07,0,1,NA
2017-09-08,1,0,NA
2017-09-09,1,0,NA
2017-09-10,1,1,NA
2017-09-11,0,0,NA
2017-09-12,0,1,NA
2017-09-13,1,1,NA
2017-09-14,1,0,NA
2017-09-15,0,1,NA
2017-09-16,0,0,NA
2017-09-17,0,1,NA
2017-09-18,1,1,NA
2017-09-19,1,0,NA
2017-09-20,0,0,NA
2017-09-21,1,1,NA
2017-09-22,1,0,NA
2017-09-23,0,1,NA
2017-09-24,0,0,NA
2017-09-25,1,1,NA
2017-09-26,0,0,NA
2017-09-27,1,0,NA
2017-09-28,0,1,NA
2017-09-29,1,1,NA
2017-09-30,0,0,NA
2017-10-01,1,0,NA
2017-10-02,0,1,NA
2017-10-03,0,0,NA
2017-10-04,1,1,NA
2017-10-05,0,1,NA
2017-10-06,1,0,NA
2017-10-07,1,1,NA
2017-10-08,1,0,NA
2017-10-09,0,0,NA
2017-10-10,0,0,NA
2017-10-11,0,1,NA
2017-10-12,0,0,NA
2017-10-13,0,1,NA
2017-10-14,0,1,NA
2017-10-15,0,1,NA
2017-10-16,0,1,NA
2017-10-17,0,1,NA
2017-10-18,0,1,NA
2017-10-19,0,1,NA
2017-10-20,0,1,NA
2017-10-21,0,1,NA
2017-10-22,0,1,NA
2017-10-23,0,1,NA
2017-10-24,0,1,NA
2017-10-25,0,1,NA
2017-10-26,0,1,NA
2017-10-27,0,1,NA
2017-10-28,0,1,NA
2017-10-29,0,1,NA
2017-10-30,0,1,NA
2017-10-31,0,1,NA
2017-11-01,0,1,NA
2017-11-02,0,1,NA ## I want an interaction effect on `Fan*Door`, since they may have different effects
## in the presence of each other. lavaan/blavaan do not support interaction effects.
## One work around, like with quadratic terms, is to create an interaction variable `Fan.r.Door.r`
## in the dataset which can be used directly in the SEM.
library(semTools)
door <- indProd(data=door, var1=2, var2=3) mp <- read.csv("~/selfexperiment/mp.csv", colClasses=c("Date", "integer")) ## Mnemosyne is a measure of cognitive performance; but taking the mean 'grade' per day doesn't work well
## since harder cards will tend to come up more often and a low daily grade may simply mean some hard cards
## came up, not that one's memory/energy is worse or better. Grades are nested within each card, but also nested
## within each day as random effects. So we do a multilevel model inferred daily random effects on grade, using
## as covariates Mnemosyne's own estimate of each card's current difficulty and how long since the last review/grade
## of each card. Then the daily random effects represent how much better than predicted based on card characteristics
## I performed each day, and hopefully a measure of how well I was thinking that day.
library(sqldf)
target <- "~/.local/share/mnemosyne/default.db"
mnemosyneGrades <- sqldf("SELECT timestamp,object_id,grade,easiness,actual_interval FROM log WHERE event_type==9;", dbname=target, method = c("integer", "factor","integer","numeric","integer"))
mnemosyneGrades$timestamp <- as.POSIXct(mnemosyneGrades$timestamp, origin = "1970-01-01", tz = "EST")
colnames(mnemosyneGrades) <- c("Timestamp", "ID", "Grade", "Easiness", "Interval.length")
mnemosyneGrades$Interval.length.log <- log1p(mnemosyneGrades$Interval.length)
mnemosyneGrades$Date <- as.Date(mnemosyneGrades$Timestamp)
library(blme)
b <- blmer(Grade ~ (1|Date) + (1|ID) + Easiness + Interval.length.log, data=mnemosyneGrades); summary(b)
dailyscores <- ranef(b)$Date
mnemosynePerformance <- data.frame(Date=as.Date(rownames(dailyscores)), Mnemosyne.score= dailyscores$"(Intercept)") night <- merge(mnemosynePerformance, merge(mp, merge(door, merge(zeo, airNight, all=TRUE), all=TRUE), all=TRUE), all=TRUE)
night[night$Date<as.Date("2016-06-02") & is.na(night$Door.r),]$Door.r <- 0 # I have always slept with door shut/fan off, so impute any missing data to 0/shut
night[night$Date<as.Date("2016-06-02") & is.na(night$Fan.r),]$Fan.r <- 0 night[night$Date<as.Date("2016-11-17") & is.na(night$Mattress.topper),]$Mattress.topper <- 0
night[night$Date>=as.Date("2016-11-17") & is.na(night$Mattress.topper),]$Mattress.topper <- 1 night[night$Date<as.Date("2016-10-13") & is.na(night$Pillow.new),]$Pillow.new <- 0
night[night$Date>=as.Date("2016-10-13") & is.na(night$Pillow.new),]$Pillow.new <- 1 night$Total.Z.2 <- night$Total.Z^2
night$CO2.N.mean.root <- sqrt(night$CO2.N.mean)
night$Time.in.Wake.log <- log1p(night$Time.in.Wake)
night$ZQ.2 <- night$ZQ^2
night$Time.in.REM.2 <- night$Time.in.REM^2
night$Time.in.Light.2 <- night$Time.in.Light^2 ## plot all variables to look for skew:
# library(reshape2)
# library(ggplot2)
# d <- melt(night[-1])
# ggplot(d,aes(x = value)) + facet_wrap(~variable,scales = "free") + geom_histogram() # write.csv(night, file="~/wiki/docs/zeo/2016-11-02-co2.csv", row.names=FALSE)
# night <- read.csv("~/wiki/docs/zeo/2016-11-02-co2.csv", colClasses=c("Date", rep("numeric", 26))) Example for single-indicator measurement error: x <- rnorm(500000)
x_measured <- x + rnorm(500000, sd=0.2)
cor(x, x_measured)
# [1] 0.9801406332
cov(x, x_measured)
# [1] 0.9942188471
## standardized, so residual variance is just 1-0.98 or ~0.02
y <- x_measured + rnorm(50000, sd=0.5)
library(lavaan)
model <- 'y ~ X
 X =~ x
 x ~~ 0.0207738304*x'
summary(sem(model, data=data.frame(x=x_measured, y=y))) library(blavaan)
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  SLEEP =~ ZQ.2_latent + Total.Z.2_latent + Time.in.REM.2_latent + Time.in.Light.2_latent + Time.in.Deep_latent + Time.to.Z.log_latent + Time.in.Wake.log_latent + Awakenings_latent + Morning.Feel  HUMIDITY =~ Humidity.K.mean + Humidity.K.max + Humidity.N.mean + Humidity.N.max
 TEMPERATURE =~ Temperature.K.mean + Temperature.K.max + Temperature.N.mean + Temperature.N.max
 CO2 =~ CO2.N.mean.root + CO2.N.max
 NOISE =~ Noise.N.mean + Noise.N.max  CO2 ~ Fan.r + Door.r + Fan.r.Door.r
 HUMIDITY ~ Fan.r + Door.r + Fan.r.Door.r
 TEMPERATURE ~ Fan.r + Door.r + Fan.r.Door.r  SLEEP ~ CO2 + HUMIDITY + TEMPERATURE + NOISE + CO2 + Fan.r + Door.r + Fan.r.Door.r + Pillow.new + Mattress.topper + Mattress.r
 MP ~ CO2 + SLEEP + Pillow.new + Mattress.topper + Mattress.r
 Mnemosyne.score ~ CO2 + SLEEP
 '
# library(lavaan)
# s1 <- sem(model1, missing="FIML", data=scale(night[-1])); summary(s1) system.time(s1 <- bsem(model1, test="none", n.chains=8, burnin=20000, sample=6000, dp = dpriors(nu = "dnorm(0,1)", alpha = "dnorm(0,1)", beta = "dnorm(0,200)"), jagcontrol=list(method="rjparallel"), fixed.x=FALSE, data=scale(night[-1]))); summary(s1)