load packages

packages <- c('tidyverse',     # data handling
              'lmerTest', 
              'modelr',
              'sjPlot',
              'kableExtra',
              'broom.mixed',
              'gridExtra',
              'brms',
              'bayesplot',
              'rstanarm',
              'performance',
              'tidybayes',
              'brant')    
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages], repos = c(CRAN = "https://cran.rstudio.com"))
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))

packv <- NULL
for (i in 1:length(packages)) {
  packv = rbind(packv, c(packages[i], as.character(packageVersion(packages[i]))))
}
colnames(packv) <- c("Package", "Version") 
packv %>% 
  as_tibble %>% 
  arrange(Package) %>% 
  kable(align = "l") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Package	Version
bayesplot	1.9.0
brant	0.3.0
brms	2.16.3
broom.mixed	0.2.9.3
gridExtra	2.3
kableExtra	1.3.4
lmerTest	3.1.3
modelr	0.1.8
performance	0.9.0
rstanarm	2.21.1
sjPlot	2.8.10
tidybayes	3.0.2
tidyverse	1.3.1

read data

# according to subject notes
# vret codes for treatment: 1= vret, 0 = ivet
file <- 'data/VR_spider_EEG-EG_final_with_clinical_data.tsv'
listEGnotes <- read.csv(file, sep = '\t', header = T)
listEGnotes$sess1 <- 1
listEGnotes$sess1[listEGnotes$pre=='nope'] <- 0
listEGnotes$sess2 <- 1
listEGnotes$sess2[listEGnotes$post=='nope'] <- 0
df <- listEGnotes
rm(listEGnotes)

prepare data

df <- df %>%
  pivot_longer(cols = BAT_before:FSQ_after, names_to = 'measure', values_to = 'rating') %>%  
  mutate(treat = ifelse(vret == 1, "VRET", "IVET"),
         treat = factor(treat, levels = c('IVET', 'VRET'))) %>% 
  mutate(sess = ifelse(str_detect(measure, "before"), "pre", "post"),
         sess = factor(sess, levels = c('pre', 'post')),
         measure = str_remove_all(measure, "_before"),
         measure = str_remove_all(measure, "_after")) %>% 
  select(-c(sess1, sess2, pre, post, vret)) %>% 
  relocate(treat, sess, measure, fp, rating) 

# define contrasts
contrasts(df$treat) <- c(-0.5, +0.5)
contrasts(df$sess) <- c(0, 1)

show data

df %>% 
  tibble()

means

tables show means (SD in parentheses).

by IVET and VRET

df %>% 
  group_by(treat, sess, measure) %>%
  summarise(N = n(),
            m_sd = sprintf("%.2f (%.2f)", 
                           mean(rating, na.rm = TRUE), 
                           sd(rating, na.rm = TRUE)), 
                           .groups = 'drop') %>% 
  # N is slightly misleading because some subjects are excluded (but still counted in N)
  pivot_wider(names_from = measure, values_from = m_sd) %>%
  rename("Treatment" = treat,
         "Session" = sess) %>%
  arrange(Treatment, Session) %>%
  kable() %>%
  kable_styling()

Treatment	Session	N	BAT	FSQ	SPQ
IVET	pre	33	5.91 (2.23)	97.03 (12.40)	22.27 (2.54)
IVET	post	33	10.84 (1.59)	50.97 (16.49)	12.64 (5.18)
VRET	pre	37	4.62 (2.84)	96.81 (16.74)	22.54 (4.23)
VRET	post	37	8.58 (2.48)	67.97 (22.93)	16.19 (6.11)

across IVET and VRET

df %>% 
  group_by(sess, measure) %>%
  summarise(N = n(),
            m_sd = sprintf("%.2f (%.2f)", 
                           mean(rating, na.rm = TRUE), 
                           sd(rating, na.rm = TRUE)), 
                           .groups = 'drop') %>% 
  # N is slightly misleading because some subjects are excluded (but still counted in N)
  pivot_wider(names_from = measure, values_from = m_sd) %>%
  rename("Session" = sess) %>%
  arrange(Session) %>%
  kable() %>%
  kable_styling()

Session	N	BAT	FSQ	SPQ
pre	70	5.23 (2.63)	96.91 (14.75)	22.41 (3.52)
post	70	9.63 (2.39)	59.84 (21.72)	14.49 (5.92)

factor definition

df %>% 
  select(where(is.factor)) %>% 
  purrr::map(contrasts)

## $treat
##      [,1]
## IVET -0.5
## VRET  0.5
## 
## $sess
##      [,1]
## pre     0
## post    1

interpretation help

intercept: average rating in the pre session across IVET and VRET

treat1: difference of VRET minus IVET

sess1: difference of post minus pre

treat1:sess1: difference of VRET minus IVET in the session effect (post minus pre)

BAT

Frequentist model

summary

data <- df %>%
  filter(measure == "BAT",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp)) 
model_name <- "clinical_data_BAT.model"

model_formula <- formula("rating ~ 1 + treat*sess + (1 | fp)")

if (file.exists(sprintf("results/models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("results/models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), lmerTest::lmer(model_formula,
                                           control = lmerControl(optimizer = "nloptwrap",
                                                                 optCtrl = list(maxfun = 1e7),
                                                                 calc.derivs = FALSE), 
                                           data = data))
  saveRDS(eval(parse(text = model_name)), sprintf("results/models/model_%s.RDS", model_name))
}
summary(get(model_name))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: model_formula
##    Data: data
## Control: 
## lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000),  
##     calc.derivs = FALSE)
## 
## REML criterion at convergence: 599
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.63355 -0.47583  0.05055  0.52444  1.72033 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  fp       (Intercept) 2.871    1.694   
##  Residual             2.672    1.635   
## Number of obs: 137, groups:  fp, 70
## 
## Fixed effects:
##              Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)    5.2654     0.2819 106.0010  18.680 <0.0000000000000002 ***
## treat1        -1.2875     0.5637 106.0010  -2.284              0.0244 *  
## sess1          4.4177     0.2817  67.0380  15.683 <0.0000000000000002 ***
## treat1:sess1  -0.9298     0.5634  67.0380  -1.650              0.1035    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) treat1 sess1 
## treat1      -0.057              
## sess1       -0.482  0.028       
## treat1:sss1  0.028 -0.482 -0.071

coef table

get(model_name) %>%
  broom.mixed::tidy(conf.int = TRUE) %>% 
  filter(effect == "fixed") %>% 
  select(-group) %>% 
  kable() %>%
  kable_styling()

effect	term	estimate	std.error	statistic	df	p.value	conf.low	conf.high
fixed	(Intercept)	5.2653563	0.2818731	18.679881	106.00095	0.0000000	4.706515	5.8241971
fixed	treat1	-1.2874693	0.5637462	-2.283774	106.00095	0.0243784	-2.405151	-0.1697876
fixed	sess1	4.4176515	0.2816871	15.682831	67.03797	0.0000000	3.855408	4.9798955
fixed	treat1:sess1	-0.9297659	0.5633742	-1.650352	67.03797	0.1035489	-2.054254	0.1947221

Bayesian model

summary

data <- df %>%
  filter(measure == "BAT",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp))
mod_brm_BAT <- brm(rating ~ 1 + treat*sess + (1 | fp),
  family = gaussian(),
  data = data,
  prior = c(prior(normal(0, 4), class = Intercept),
            prior(normal(0, 4), class = b)),
 chains = 4,
  file = "results/models/clinical_data_BAT_brm.model",
  # file to save/reuse model
  cores = 4,
  iter = 3000,
  warmup = 1000,
  seed = sum(as.integer(iconv("BAT", toRaw = TRUE)[[1]])),
  save_pars = save_pars(all = TRUE),
  sample_prior = "yes")
# Model Summary.
summary(mod_brm_BAT)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rating ~ 1 + treat * sess + (1 | fp) 
##    Data: data (Number of observations: 137) 
##   Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~fp (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.69      0.24     1.23     2.19 1.00     2275     3080
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        5.25      0.29     4.69     5.82 1.00     6632     5915
## treat1          -1.29      0.56    -2.38    -0.19 1.00     6684     6328
## sess1            4.40      0.29     3.82     4.97 1.00    17311     5422
## treat1:sess1    -0.92      0.57    -2.03     0.20 1.00    15497     6249
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.67      0.15     1.41     2.00 1.00     2939     4414
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

coef table

posterior_mod <- as.array(mod_brm_BAT$fit)
coeffsall <- dimnames(posterior_mod)$parameters %>% enframe() %>% 
  mutate(FE = str_detect(value, pattern = "b_") | value == "sigma")
tidyMCMC(mod_brm_BAT$fit,
      conf.int = TRUE,
      rhat = TRUE,
      ess = TRUE,
      conf.level = 0.95,
      conf.method = "quantile",
      pars = coeffsall$value[coeffsall$FE]) %>% 
      kable(digits = 4) %>%
      kable_styling()

term	estimate	std.error	conf.low	conf.high	rhat	ess
b_Intercept	5.2497	0.2896	4.6876	5.8176	0.9999	6594
b_treat1	-1.2850	0.5615	-2.3808	-0.1875	0.9999	6644
b_sess1	4.3966	0.2932	3.8187	4.9684	0.9999	17443
b_treat1:sess1	-0.9161	0.5672	-2.0302	0.2050	0.9997	15876
sigma	1.6595	0.1503	1.4053	1.9970	1.0013	2922

FSQ

Frequentist model

summary

data <- df %>%
  filter(measure == "FSQ",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp)) 
model_name <- "clinical_data_FSQ.model"

model_formula <- formula("rating ~ 1 + treat*sess + (1 | fp)")

if (file.exists(sprintf("results/models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("results/models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), lmerTest::lmer(model_formula,
                                           control = lmerControl(optimizer = "nloptwrap",
                                                                 optCtrl = list(maxfun = 1e7),
                                                                 calc.derivs = FALSE), 
                                           data = data))
  saveRDS(eval(parse(text = model_name)), sprintf("results/models/model_%s.RDS", model_name))
}
summary(get(model_name))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: model_formula
##    Data: data
## Control: 
## lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000),  
##     calc.derivs = FALSE)
## 
## REML criterion at convergence: 1164.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.40738 -0.58708  0.02822  0.57483  1.82766 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  fp       (Intercept) 106.4    10.32   
##  Residual             205.3    14.33   
## Number of obs: 139, groups:  fp, 70
## 
## Fixed effects:
##              Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)   96.9206     2.1135 121.1524  45.857 < 0.0000000000000002 ***
## treat1        -0.2195     4.2270 121.1524  -0.052             0.958673    
## sess1        -37.4202     2.4362  67.8811 -15.360 < 0.0000000000000002 ***
## treat1:sess1  17.2807     4.8725  67.8811   3.547             0.000713 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) treat1 sess1 
## treat1      -0.057              
## sess1       -0.571  0.033       
## treat1:sss1  0.033 -0.571 -0.048

coef table

get(model_name) %>%
  broom.mixed::tidy(conf.int = TRUE) %>% 
  filter(effect == "fixed") %>% 
  select(-group) %>% 
  kable() %>%
  kable_styling()

effect	term	estimate	std.error	statistic	df	p.value	conf.low	conf.high
fixed	(Intercept)	96.9205569	2.113516	45.8574960	121.15244	0.0000000	92.736347	101.104766
fixed	treat1	-0.2194922	4.227032	-0.0519258	121.15244	0.9586733	-8.587911	8.148927
fixed	sess1	-37.4202473	2.436234	-15.3598768	67.88111	0.0000000	-42.281830	-32.558665
fixed	treat1:sess1	17.2807176	4.872467	3.5466053	67.88111	0.0007126	7.557552	27.003883

Bayesian model

summary

data <- df %>%
  filter(measure == "FSQ",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp))
mod_brm_FSQ <- brm(rating ~ 1 + treat*sess + (1 | fp),
  family = gaussian(),
  data = data,
  prior = c(prior(normal(0, 4), class = Intercept),
            prior(normal(0, 4), class = b)),
 chains = 4,
  file = "results/models/clinical_data_FSQ_brm.model",
  # file to save/reuse model
  cores = 4,
  iter = 3000,
  warmup = 1000,
  seed = sum(as.integer(iconv("FSQ", toRaw = TRUE)[[1]])),
  save_pars = save_pars(all = TRUE),
  sample_prior = "yes")
# Model Summary.
summary(mod_brm_FSQ)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rating ~ 1 + treat * sess + (1 | fp) 
##    Data: data (Number of observations: 139) 
##   Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~fp (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)    61.17      6.97    48.37    75.31 1.00     1202     2175
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       29.33      4.89    19.73    39.03 1.00     5113     5393
## treat1           0.54      3.89    -7.04     8.05 1.00     4863     4161
## sess1          -23.30      3.14   -29.17   -16.83 1.00     3577     4371
## treat1:sess1     5.00      3.42    -1.71    11.66 1.00     7662     6037
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    18.21      2.08    14.63    22.73 1.00     1820     3462
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

coef table

posterior_mod <- as.array(mod_brm_FSQ$fit)
coeffsall <- dimnames(posterior_mod)$parameters %>% enframe() %>% 
  mutate(FE = str_detect(value, pattern = "b_") | value == "sigma")
tidyMCMC(mod_brm_FSQ$fit,
      conf.int = TRUE,
      rhat = TRUE,
      ess = TRUE,
      conf.level = 0.95,
      conf.method = "quantile",
      pars = coeffsall$value[coeffsall$FE]) %>% 
      kable(digits = 4) %>%
      kable_styling()

term	estimate	std.error	conf.low	conf.high	rhat	ess
b_Intercept	29.3456	4.8873	19.7270	39.0347	1.0001	5068
b_treat1	0.5489	3.8860	-7.0387	8.0498	1.0000	4847
b_sess1	-23.4233	3.1407	-29.1743	-16.8348	1.0000	3569
b_treat1:sess1	5.0195	3.4185	-1.7141	11.6600	0.9999	7578
sigma	18.0516	2.0804	14.6288	22.7329	1.0002	1834

SPQ

Frequentist model

summary

data <- df %>%
  filter(measure == "SPQ",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp)) 
model_name <- "clinical_data_SPQ.model"

model_formula <- formula("rating ~ 1 + treat*sess + (1 | fp)")

if (file.exists(sprintf("results/models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("results/models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), lmerTest::lmer(model_formula,
                                           control = lmerControl(optimizer = "nloptwrap",
                                                                 optCtrl = list(maxfun = 1e7),
                                                                 calc.derivs = FALSE), 
                                           data = data))
  saveRDS(eval(parse(text = model_name)), sprintf("results/models/model_%s.RDS", model_name))
}
summary(get(model_name))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: model_formula
##    Data: data
## Control: 
## lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000),  
##     calc.derivs = FALSE)
## 
## REML criterion at convergence: 807.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12214 -0.50349 -0.04204  0.66607  1.78063 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  fp       (Intercept)  7.772   2.788   
##  Residual             14.542   3.813   
## Number of obs: 139, groups:  fp, 70
## 
## Fixed effects:
##              Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)   22.4066     0.5655 120.6680  39.621 <0.0000000000000002 ***
## treat1         0.2678     1.1310 120.6680   0.237              0.8132    
## sess1         -7.9938     0.6485  67.9240 -12.327 <0.0000000000000002 ***
## treat1:sess1   3.2850     1.2969  67.9240   2.533              0.0136 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) treat1 sess1 
## treat1      -0.057              
## sess1       -0.568  0.032       
## treat1:sss1  0.032 -0.568 -0.048

coef table

get(model_name) %>%
  broom.mixed::tidy(conf.int = TRUE) %>% 
  filter(effect == "fixed") %>% 
  select(-group) %>% 
  kable() %>%
  kable_styling()

effect	term	estimate	std.error	statistic	df	p.value	conf.low	conf.high
fixed	(Intercept)	22.4066339	0.5655178	39.6214497	120.66796	0.0000000	21.2870112	23.526257
fixed	treat1	0.2678133	1.1310355	0.2367859	120.66796	0.8132244	-1.9714321	2.507059
fixed	sess1	-7.9938447	0.6484745	-12.3271531	67.92398	0.0000000	-9.2878817	-6.699808
fixed	treat1:sess1	3.2850378	1.2969490	2.5328966	67.92398	0.0136274	0.6969638	5.873112

Bayesian model

summary

data <- df %>%
  filter(measure == "SPQ",
         !is.na(rating)) %>% 
  mutate(fp = as.factor(fp))
mod_brm_SPQ <- brm(rating ~ 1 + treat*sess + (1 | fp),
  family = gaussian(),
  data = data,
  prior = c(prior(normal(0, 4), class = Intercept),
            prior(normal(0, 4), class = b)),
 chains = 4,
  file = "results/models/clinical_data_SPQ_brm.model",
  # file to save/reuse model
  cores = 4,
  iter = 3000,
  warmup = 1000,
  seed = sum(as.integer(iconv("SPQ", toRaw = TRUE)[[1]])),
  save_pars = save_pars(all = TRUE),
  sample_prior = "yes")
# Model Summary.
summary(mod_brm_SPQ)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rating ~ 1 + treat * sess + (1 | fp) 
##    Data: data (Number of observations: 139) 
##   Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~fp (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.70      0.61     1.34     3.79 1.00     1725     1875
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       22.04      0.57    20.89    23.14 1.00     7217     6978
## treat1           0.43      1.08    -1.64     2.57 1.00     6218     5293
## sess1           -7.76      0.65    -9.04    -6.47 1.00    14361     6179
## treat1:sess1     2.94      1.24     0.51     5.37 1.00     9742     6304
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     3.91      0.36     3.29     4.69 1.00     2232     2653
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

coef table

posterior_mod <- as.array(mod_brm_SPQ$fit)
coeffsall <- dimnames(posterior_mod)$parameters %>% enframe() %>% 
  mutate(FE = str_detect(value, pattern = "b_") | value == "sigma")
tidyMCMC(mod_brm_SPQ$fit,
      conf.int = TRUE,
      rhat = TRUE,
      ess = TRUE,
      conf.level = 0.95,
      conf.method = "quantile",
      pars = coeffsall$value[coeffsall$FE]) %>% 
      kable(digits = 4) %>%
      kable_styling()

term	estimate	std.error	conf.low	conf.high	rhat	ess
b_Intercept	22.0480	0.5686	20.8934	23.1357	0.9996	7146
b_treat1	0.4083	1.0805	-1.6446	2.5663	1.0010	6173
b_sess1	-7.7742	0.6537	-9.0369	-6.4659	0.9997	14131
b_treat1:sess1	2.9376	1.2414	0.5146	5.3665	1.0005	9648
sigma	3.8889	0.3592	3.2895	4.6947	1.0019	2132

VR spiders: analyses of clinical data

Stefan Wiens

2022-05-10

load packages

read data

prepare data

show data

means

factor definition

interpretation help

BAT

Frequentist model

summary

coef table

Bayesian model

summary

coef table

FSQ

Frequentist model

summary

coef table

Bayesian model

summary

coef table

SPQ

Frequentist model

summary

coef table

Bayesian model

summary

coef table