Iris - VR presentation
Russell Chan
2024-06-20
# Packages
library(lme4)## Loading required package: Matrixlibrary(readxl)
library(haven)
library(effects)## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.library(lattice)
library(car)
library(ggplot2)
library(knitr)
library(reshape2)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forcats)
library(DHARMa)## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(phia)
library(lsmeans)## Loading required package: emmeans## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.library(emmeans)
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:Hmisc':
##
## subplot## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(tinytex)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ✔ purrr 1.0.2 ✔ tibble 3.2.1
## ✔ readr 2.1.5 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ plotly::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ plotly::select() masks MASS::select(), dplyr::select()
## ✖ purrr::some() masks car::some()
## ✖ Hmisc::src() masks dplyr::src()
## ✖ Hmisc::summarize() masks dplyr::summarize()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(devtools)## Loading required package: usethis
##
## Attaching package: 'devtools'
##
## The following object is masked from 'package:emmeans':
##
## testlibrary(brms)## Loading required package: Rcpp
## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
##
## The following object is masked from 'package:survival':
##
## kidney
##
## The following object is masked from 'package:lme4':
##
## ngrps
##
## The following object is masked from 'package:stats':
##
## arlibrary(ggthemes)
library(readxl)
library(janitor)##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractlibrary(tidyverse)
library(mascutils)##
## Attaching package: 'mascutils'
##
## The following object is masked from 'package:tidyr':
##
## expand_grid
##
## The following object is masked from 'package:car':
##
## logit
##
## The following object is masked from 'package:base':
##
## modelibrary(bayr)## Registered S3 methods overwritten by 'bayr':
## method from
## coef.brmsfit brms
## knit_print.tbl_obs mascutils
## predict.brmsfit brms
## print.tbl_obs mascutils
##
## Attaching package: 'bayr'
##
## The following objects are masked from 'package:mascutils':
##
## as_tbl_obs, discard_all_na, discard_redundant, expand_grid,
## go_arrange, go_first, left_union, reorder_levels, rescale_centered,
## rescale_unit, rescale_zero_one, update_by, z_score, z_trans
##
## The following objects are masked from 'package:brms':
##
## fixef, ranef
##
## The following object is masked from 'package:tidyr':
##
## expand_grid
##
## The following objects are masked from 'package:lme4':
##
## fixef, raneflibrary(brms)
library(ggthemes)
library(afex)## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
##
## The following object is masked from 'package:lme4':
##
## lmer#SetWD
setwd("/Users/rc/Documents/Iris_EDA")library(readxl)
PhysioDF <- read_excel("physiological-data-merged.xlsx", col_types = c("numeric", "numeric", "numeric", "numeric", "numeric", "numeric"))PhysioDF$Subject = factor(PhysioDF$`Subject no.`)
PhysioDF$Session = factor(PhysioDF$Session, ordered = TRUE, levels=c('1' ,'2', '3'))
PhysioDF$Group = factor(PhysioDF$Group)
PhysioDF$Minute = factor(PhysioDF$Minute)
PhysioDF$BPM = as.numeric(PhysioDF$`BPM(PR)`)BPM <- lmer(BPM ~ Group * Session * Minute + (1|Subject), data = PhysioDF, REML = FALSE)
Anova(BPM)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: BPM
## Chisq Df Pr(>Chisq)
## Group 0.0816 1 0.7752
## Session 3.6320 2 0.1627
## Minute 21.6062 2 2.034e-05 ***
## Group:Session 3.0554 2 0.2170
## Group:Minute 0.3612 2 0.8348
## Session:Minute 4.9910 4 0.2882
## Group:Session:Minute 1.0732 4 0.8985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(BPM, ddf="Satterthwaite")## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: BPM ~ Group * Session * Minute + (1 | Subject)
## Data: PhysioDF
##
## AIC BIC logLik deviance df.resid
## 3424.7 3506.2 -1692.4 3384.7 414
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5222 -0.4719 0.0414 0.4542 4.7465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 98.24 9.911
## Residual 111.69 10.569
## Number of obs: 434, groups: Subject, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 83.3660 2.3382 73.9510 35.654 < 2e-16 ***
## Group2 0.8562 3.3334 73.3206 0.257 0.79800
## Session.L -0.6915 2.1648 385.6795 -0.319 0.74956
## Session.Q -2.4422 2.1309 385.2604 -1.146 0.25247
## Minute2 4.7011 1.7521 385.2226 2.683 0.00761 **
## Minute3 5.3868 1.7449 385.1276 3.087 0.00217 **
## Group2:Session.L 3.9030 3.0562 385.3613 1.277 0.20234
## Group2:Session.Q 0.6221 3.0322 385.1490 0.205 0.83757
## Group2:Minute2 0.7059 2.4896 385.1846 0.284 0.77691
## Group2:Minute3 -0.7793 2.4845 385.1376 -0.314 0.75393
## Session.L:Minute2 -2.5621 3.0573 385.3097 -0.838 0.40255
## Session.Q:Minute2 0.8587 3.0121 385.1328 0.285 0.77574
## Session.L:Minute3 -2.6025 3.0387 385.1698 -0.856 0.39229
## Session.Q:Minute3 0.3268 3.0058 385.0844 0.109 0.91348
## Group2:Session.L:Minute2 -1.3655 4.3323 385.2548 -0.315 0.75279
## Group2:Session.Q:Minute2 1.0687 4.2917 385.1131 0.249 0.80348
## Group2:Session.L:Minute3 -3.2901 4.3060 385.1047 -0.764 0.44530
## Group2:Session.Q:Minute3 2.9821 4.3005 385.1705 0.693 0.48846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itEDA <- lmer(EDA ~ Group * Session * Minute + (1|Subject), data = PhysioDF, REML = FALSE)
Anova(EDA)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: EDA
## Chisq Df Pr(>Chisq)
## Group 1.2304 1 0.267320
## Session 10.8460 2 0.004414 **
## Minute 4.4329 2 0.108997
## Group:Session 1.1983 2 0.549283
## Group:Minute 0.8803 2 0.643943
## Session:Minute 4.4645 4 0.346782
## Group:Session:Minute 2.1819 4 0.702339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(EDA, ddf="Satterthwaite")## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: EDA ~ Group * Session * Minute + (1 | Subject)
## Data: PhysioDF
##
## AIC BIC logLik deviance df.resid
## 681.5 763.3 -320.8 641.5 421
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6331 -0.2408 0.0159 0.1509 12.1724
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.2214 0.4705
## Residual 0.1913 0.4374
## Number of obs: 441, groups: Subject, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.26640 0.10680 67.42035 2.494 0.0151 *
## Group2 -0.11071 0.15260 67.42035 -0.725 0.4707
## Session.L 0.08853 0.08748 392.00000 1.012 0.3122
## Session.Q -0.19841 0.08748 392.00000 -2.268 0.0239 *
## Minute2 0.05907 0.07143 392.00000 0.827 0.4088
## Minute3 0.15293 0.07143 392.00000 2.141 0.0329 *
## Group2:Session.L -0.03668 0.12500 392.00000 -0.293 0.7694
## Group2:Session.Q 0.09260 0.12500 392.00000 0.741 0.4592
## Group2:Minute2 -0.04087 0.10206 392.00000 -0.400 0.6890
## Group2:Minute3 -0.09543 0.10206 392.00000 -0.935 0.3503
## Session.L:Minute2 -0.13294 0.12372 392.00000 -1.074 0.2833
## Session.Q:Minute2 0.08263 0.12372 392.00000 0.668 0.5046
## Session.L:Minute3 -0.28284 0.12372 392.00000 -2.286 0.0228 *
## Session.Q:Minute3 0.13195 0.12372 392.00000 1.066 0.2869
## Group2:Session.L:Minute2 0.11820 0.17678 392.00000 0.669 0.5041
## Group2:Session.Q:Minute2 -0.09147 0.17678 392.00000 -0.517 0.6051
## Group2:Session.L:Minute3 0.22274 0.17678 392.00000 1.260 0.2084
## Group2:Session.Q:Minute3 -0.13297 0.17678 392.00000 -0.752 0.4524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# plot BPM
ae.m.BPM <- allEffects(BPM)
ae.m.df.BPM = as.data.frame(ae.m.BPM[[1]])
ae.BPM = ggplot(ae.m.df.BPM,
aes(x = Minute, y = fit, group = Group, color = Group)) +
geom_errorbar(aes(ymin = fit-se, ymax = fit+se), width=.1) +
geom_line() +
geom_point() +
ylab("Heart Rate (BPM)") +
xlab("Minute") +
facet_grid(~Session)+
theme_classic()
plot(ae.BPM)
# plot BPM
ae.m.EDA <- allEffects(EDA)
ae.m.df.EDA = as.data.frame(ae.m.EDA[[1]])
ae.EDA = ggplot(ae.m.df.EDA,
aes(x = Minute, y = fit, group = Group, color = Group)) +
geom_errorbar(aes(ymin = fit-se, ymax = fit+se), width=.1) +
geom_line() +
geom_point() +
ylab("Skin Conductance (EDA)") +
xlab("Minute") +
facet_grid(~Session)+
theme_classic()
plot(ae.EDA)
aov_BPM <- aov(BPM ~ Group * Session * Minute + Error(Subject), data = PhysioDF)
summary(aov_BPM)##
## Error: Subject
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 115 114.9 0.106 0.746
## Session 2 1438 718.8 0.664 0.520
## Minute 2 275 137.3 0.127 0.881
## Residuals 43 46546 1082.5
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Session 2 408 204.1 1.751 0.175
## Minute 2 2412 1205.9 10.347 4.24e-05 ***
## Group:Session 2 350 175.0 1.502 0.224
## Group:Minute 2 39 19.6 0.168 0.845
## Session:Minute 4 553 138.3 1.186 0.316
## Group:Session:Minute 4 122 30.5 0.262 0.902
## Residuals 369 43004 116.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov_EDA <- aov(EDA ~ Group * Session * Minute + Error(Subject), data = PhysioDF)
summary(aov_EDA)##
## Error: Subject
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 2.69 2.687 1.18 0.283
## Residuals 47 106.99 2.276
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Session 2 2.08 1.0376 5.202 0.00591 **
## Minute 2 0.85 0.4241 2.126 0.12075
## Group:Session 2 0.23 0.1146 0.575 0.56337
## Group:Minute 2 0.17 0.0842 0.422 0.65593
## Session:Minute 4 0.85 0.2135 1.071 0.37078
## Group:Session:Minute 4 0.42 0.1044 0.523 0.71873
## Residuals 376 75.00 0.1995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1