QS Experiments / QS - Which factors influence my sleep

Abstract

What did i do?

I've measured my food intake, physical activity, positive and negative affects, sickness, sleep quality and quantity to find out how different factors affect my sleep.

How did i do it?

I've measured all variables of interest for ~14 months and there is 115 days with all data present. Days without all metrics present were excluded.

What did i learn?

Caffeine, bedtime start, sickness, vitamin D3 and negative emotions statistically significantly and independently influence my sleep. I'll go into direction of lowering caffeine intake, going sleep earlier, take D3 every other morning, and try to reduce stress. New findings going in same direction with previous findings [1] [2].

Introduction

There is a lof of lifestyle factors which may disrupt or improve sleep quantity and quality.

The purpose of this observational data analysis (n=1) is to determine factors influencing my total sleep time.

Materials & Methods

Participants

1 adult male anthropometrics was described in previous article

Experimental design

During 115 nights (from 2020-08-27 to 2021-09-16) sleep quality and quantity was assessed by Dreem 2 EEG Headband which was validated against gold standard PSG. At same time, all food and fluid intake was carefully logged and measures of intereset were extracted from dietary log. Content of nutrients in each product was determined from a popular online food database - https://nutritiondata.self.com/
Measures of interest: daily calories, protein (cal), fat (cal), carb (cal), animal products (meat, fish, eggs in g), d3 (mcg), b12 (μg), potassium (mcg), magnesium (mcg).

Physical activity was calculated by using sRPE approach, described here.

Positive and Negative affect was measured using PANAS scale, described here.

Sauna and occasional meditation sessions were measured in minutes.

Subjective sleep quality, mood, depression, anxiety, stress, fatigue and sleepiness was assessed by a popular in scientific research questionnaires described in previous experiment

Results

Full dataset summary:

Running simple linear regression with all params versus total sleep time (tst):

After checking, i've decided to remove non-significant variables. Since d3 isnt far from significance (p-value 0.127), i've decided to keep it. Physical activity, food nutrient structure, b12, potassium & magnesium didnt reveal any connection with a sleep time. That might be because most of my activity falls on the morning / first half of day and i dont eat 3-5 hours prior bed which results in not detectable (small) effects by rough linear analysis or being masked by other influencers not taken into account.

Linear regression after removing non-significant variables:

We can see improvements in model: F-statistic increased and R-squared also increased, explaining 23.2% of total sleep time variability.

Linear regression results, Total Sleep Time

effect, minutes p-adjusted
Negative Affect -5.090 0.009* Worse
Sick score 8.190 0.00012* Better
Vitamin D3, per 100 mcg 6.121 0.0176* Better
Bedtime start, h from 00:00 -15.498 0.0424* Worse
Caffeine, per 100 mg -29.552  0.0017* Worse

Finally, this model tells:

Discussion

The main result of this experiment is a statistically significant and independent associations between sleep in restorative phases and caffeine, bedtime start, vitamin D3, sickness, negative affect.
Results pointing me to

Checks for linearity, residuals distribution, interaction, collinearity, homoscendasicity, error term correlation seems to indicate acceptable quality of model and explained in Statystical Analysis section.

A limitations:

Statistical Analysis

P-values surviving adjusting for multiple comparisons:

The regression lines for each predictor:

95% Confidence intervals for slopes are narrow enough and Effect size for model is moderate, based on adjusted R-squared.

Residuals and their distribution (QQ plot) looks acceptable

Median of residuals is 8 minutes, not far from 0 and looks acceptable

Model passed test for heteroscendacity, confirms acceptable constant error variance :

Lets check for collinearity in variables of interest:

We can see a correlations between sick and negative affect and between d3 and starttime. Lets try to account for their potential interaction and see if model is improved or explained something new: (tst ~ na + sick + sick*na + d3 + starttime + d3*starttime + caff) results in F-statistic 6.847 and shows no improvement, indicating that (tst ~ na + sick + d3 + starttime + caff) with F = 9 is better.
Lets do a final statistical test for collinearity by checking if VIF (variance inflation factor) for each predictor is less than 10

VIF test makes me confident in abscense of collinearity which is good sign for a model.

We can see slight autocorrelation for error term, but it looks acceptable

Data availability & Information

Welcome for questions, suggestions and critics in comments below.

Dreem 2 device unmodified data is fully available here and format is explained by manufacturer.
Daily data available here.

R Code:

library(effectsize)
library(splusTimeDate)
library(lubridate)
library(ggplot2)
library(foreach)
library(iterators)
library(skimr)
library(parallel)
library(doParallel)
library(PerformanceAnalytics)
library(car)
library(mctest)
library(ppcor)
library(corpcor)
library(lmtest)

library(mise)
mise()

registerDoParallel(cores=24)
options(scipen = 999)
alpha <- 0.05

daily <- read.csv("https://blog.kto.to/uploads/sleepv2-daily.csv")
#daily <- daily[!is.na(daily$`cal`),]
#daily <- daily[!is.na(daily$`pa`),]
#daily <- daily[!is.na(daily$`rpe`),]
daily <- daily[!is.na(daily$`na`),]
daily <- daily[!is.na(daily$`sick`),]
#daily <- daily[!is.na(daily$`food`),]
#daily <- daily[!is.na(daily$`steps`),]
daily$`sauna`[is.na(daily$`sauna`)] <- 0;
daily$`animalprot`[is.na(daily$`animalprot`)] <- 0;
daily$`mag`[is.na(daily$`mag`)] <- 0;
daily$`pot`[is.na(daily$`pot`)] <- 0;
daily$`d3`[is.na(daily$`d3`)] <- 0;
daily$`caff`[is.na(daily$`caff`)] <- 0;
daily$`b12`[is.na(daily$`b12`)] <- 0;
daily$`meditation`[is.na(daily$`meditation`)] <- 0;
#row.names(daily) <- NULL

dreem <- read.csv("https://blog.kto.to/uploads/sleepv2-dreem.csv", skip = 5, sep = ';', header = TRUE)
dreem <- dreem[!is.na(dreem$Type),]

results <- foreach(i=1:nrow(dreem), .combine='rbind', .multicombine=TRUE, .packages = "lubridate") %dopar% {
  sleepstart <- strsplit(dreem$Start.Time[i], "T")[[1]][2]
  sleepstartdate <- strsplit(dreem$Start.Time[i], "T")[[1]][1]
  resulted = c();
  for(j in 1:nrow(daily))
  {
    #if((as.Date(daily$datetime[j]) +days(1)) == as.Date(sleepstartdate))
    if((as.Date(daily$datetime[j])) == as.Date(sleepstartdate))
    {
      sleepstarttime <- strsplit(sleepstart, "\\+")[[1]][1]
      if((!is.na(daily$steps[j])))
      {
        resulted <- c(resulted, c(
          starttime = period_to_seconds(hms(sleepstarttime))/60,
          deep = period_to_seconds(hms(dreem$Deep.Sleep.Duration[i]))/60,
          rem = period_to_seconds(hms(dreem$REM.Duration[i]))/60,
          light = period_to_seconds(hms(dreem$Light.Sleep.Duration[i]))/60,
          tst = (period_to_seconds(hms(dreem$Deep.Sleep.Duration[i])) + period_to_seconds(hms(dreem$REM.Duration[i])) + period_to_seconds(hms(dreem$Light.Sleep.Duration[i])))/60,
          sol = period_to_seconds(hms(dreem$Sleep.Onset.Duration[i]))/60,
          waso  = period_to_seconds(hms(dreem$Wake.After.Sleep.Onset.Duration[i]))/60,
          awakes  = dreem$Number.of.awakenings[i],
          awake = period_to_seconds(hms(dreem$Wake.After.Sleep.Onset.Duration[i]))/60 + period_to_seconds(hms(dreem$Wake.After.Sleep.Onset.Duration[i]))/60,
          stims = dreem$Number.of.Stimulations[i],
          positions = dreem$Position.Changes[i],
          eff = dreem$Sleep.efficiency[i],
          steps = daily$steps[j]/1000,
          pa = daily$pa[j],
          sauna = daily$sauna[j],
          meditation = daily$meditation[j],
          na = daily$na[j],
          food = daily$food[j]/1000,
          rpe = daily$rpe[j]/100,
          rpe7d = daily$rpe7d[j]/100,
          sick = daily$sick[j],
          zq = 8.5*((((period_to_seconds(hms(dreem$Deep.Sleep.Duration[i])) + period_to_seconds(hms(dreem$REM.Duration[i])) + period_to_seconds(hms(dreem$Light.Sleep.Duration[i])))/3600)*1) + (0.5*(period_to_seconds(hms(dreem$REM.Duration[i]))/3600) + 1.5*(period_to_seconds(hms(dreem$Deep.Sleep.Duration[i]))/3600)) - ((period_to_seconds(hms(dreem$Wake.After.Sleep.Onset.Duration[i]))/3600+period_to_seconds(hms(dreem$Sleep.Onset.Duration[i]))/3600)*0.5 + dreem$Number.of.awakenings[i]/15)),
          d3 = daily$d3[j],
          b12 = daily$b12[j],
          cal = daily$cal[j],
          caff = daily$caff[j],
          prot = daily$prot[j],
          animalprot = daily$animalprot[j],
          carb = daily$carb[j],
          fat = daily$fat[j],
          pot = daily$pot[j],
          mag = daily$mag[j],
          startdate = as.Date(sleepstartdate)
        ))
      }
    }
  }
  return(c(resulted))
}

res_li <- as.data.frame(results)
row.names(res_li) <- NULL
summary(res_li)

res_li$d3 <- res_li$d3/100
res_li$pot <- res_li$pot/1000
res_li$caff <- res_li$caff/100
res_li$starttime <- res_li$starttime/60

#plat responce variable
plot(as_date(res_li$startdate, origin = lubridate::origin), res_li$tst)

#l <- lm(tst ~ na + sick + d3 + starttime + caff + b12 + cal + prot + carb + fat + rpe + pa + animalprot + pot + mag + sauna + meditation, data=res_li); s <- summary(l); s
l <- lm(tst ~ na + sick + d3 + starttime + caff, data=res_li); s <- summary(l); s
summary(anova(l))
mse <- mean(residuals(l)^2); mse
rss <- sum(residuals(l)^2); rss
rse <- sqrt(rss / l$df.residual); rse

#Effect size for R^2
interpret_r2(summary(l)$adj.r.squared[1], rules = "cohen1988")

#1 adjust p-values for multiple comparisons
p.adjust(summary(l)$coefficients[,4], method="BH") < alpha
confint(l, level = alpha)

print("1. Verification of Lineriaty: check each predictor")
#inspect each predictor
avPlots(l)

print("2. Verification of interaction: remove additive/interaction assumption")
#if F/T-statistics are improved, there might be an interaction
#check for interactions p-values. multiplication and both interaction predictors must be included!
la <- lm(tst ~ na + sick + (na * sick) + d3 + starttime + (d3 * starttime) + caff, data=res_li); s <- summary(la); s
la <- lm(tst ~ na + sick + d3 + starttime + (d3 * starttime) + caff, data=res_li); s <- summary(la);  s
la <- lm(tst ~ na + sick + (na * sick) + d3 + starttime + caff, data=res_li); s <- summary(la);  s

print("3. Verification of Residuals: check linearity of residuals and non-constant error variance (heteroscedasticity)")
#Linear fit. Red residual line should be linear horizontal on Residual vs Fitted
#Homoscedasticity : Residual vs Fitted and Scale-Location. Completely random, equal distribution is required to confirm Homoscedasticity
par(mfrow = c(2,2)); plot(l)
car::ncvTest(l)$p < alpha #TRUE indicates Homoscedasticity

print("4. Verification of Error term correlation")
acf(l$residuals, main="Error term correlation")
acf(l$residuals, main="Error term partial correlation")

print("5. Remove Outliers and re-check model")
cooksD <- cooks.distance(l)
influential <- cooksD[(cooksD > (20 * mean(cooksD, na.rm = TRUE)))]
names_of_influential <- names(influential)
res_li_no_outliers <- res_li[-c(as.numeric(names(influential))),]
lno <- lm(tst ~ na + sick + d3 + starttime + caff, data=res_li_no_outliers); s <- summary(lno); s
p.adjust(summary(lno)$coefficients[,4], method="BH") < alpha
par(mfrow = c(2,2)); plot(lno)
car::ncvTest(lno)$p < alpha #TRUE indicates Homoscedasticity
l <- lno

print("6. Verification of Collinearity: predictors shouldn't have collinearity")
print("Collinearity #1: visual examination (paired)")

#look at correlation chart
rl <- list(res_li$na, res_li$sick, res_li$d3, res_li$starttime, res_li$caff)
names(rl) <- c("na","sick","d3","starttime","caff")
chart.Correlation(as.data.frame(rl), histogram = TRUE, method = "pearson")

#look p-val < alpha
print("Collinearity #2 check for p-values for correlations")
pvalcor <- pcor(as.data.frame(rl), method = "pearson"); diag(pvalcor$p.value) <- NA;  diag(pvalcor$estimate) <- NA
round(pvalcor$estimate,2)
pvalcor$p.value < alpha

#look for partical correlation table
print("Collinearity #3: check partial correlation which take into account other predictors")
partialcor <- cor2pcor(cov(as.data.frame(rl))); diag(partialcor) <- NA
round(partialcor,2)

print("Collinearity #4: check variance inflation factor (for MLM). if VIF > 10 then collinearity is detected")
omcdiag(l, Inter = TRUE)
imcdiag(l, method = "VIF")$idiags[,1] > 10

Statistical analysis

RStudio version 1.3.959 and R version 4.0.2 was user for a multiple linear regression model and to calculate slopes and p-values.
P-adjusted is p-value adjusted for multiple comparisons by method of Benjamini, Hochberg, and Yekutieli.
Effect sizes based on adjusted R2, Cohen's 1988 rules
TODO: deeper check for correlation of error terms, outliers, non-constant variance of error terms, high-leverage points and collinearity.