Project.Rmd

---
title: " Appendix: Leistungsnachweis: fortgeschrittene Datenanalyse mit R"
author: "Mirko Bristle"
date: "27 10 2017"
output: word_document

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,eval=F)
```


# Setup: 
Always run this Chunk first.
```{r eval=F, include=T}
#create Operating System Path prefix
if (Sys.info()['sysname'] == "Darwin" ) {
  SysDir="/Volumes/SNFAue/"} else {
    SysDir = "N:/"}
# load necessary library files
if (!require("pacman")) install.packages("pacman")
library(pacman)
p_load(readr)
p_load(lme4)
p_load(tableone)
p_load(tidyverse)
p_load(ggplot2)
p_load(merTools)
p_load(brms)
p_load(shinystan)
p_load(ggeffects)
p_load(ggeffects)


```

### Einlesen der physiologischen Daten:

```{r eval=F, include=T}
options(max.print = 99999999)
PhysioData_org<-read_rds(paste0(SysDir,"SoOp/deprecated/Soccer deprecated/original/analysis/ECG/AggregatedData/PhysioData_org_noArt_3sdMarked_.rds"))

#exclude Bsp trails
PhysioData<-PhysioData_org  %>% filter(!grepl("Bsp",Spielfeld))
#create condition
PhysioData$Factor1<-if_else(
  grepl(paste(c("Selbst","Ingroup1","Ingroup2"),collapse="|"),PhysioData$Gruppe),
  "Us","Them")
PhysioData$Factor2<-if_else(grepl(paste(c("Selbst","Konkurrent"),
                                        collapse="|"),PhysioData$Gruppe),"Single",
                            if_else(grepl(paste(c("Ingroup1","Outgroup1"),
                                                collapse="|"),PhysioData$Gruppe),
                                    "Team_1","Team_2")) 

#make Participant ID unique by adding the between group to the name
PhysioData$BiopacSubject<-paste0(
  if_else(PhysioData$betweenCond=="soccerPlayer","sp","nsp"),
  "_",PhysioData$BiopacSubject) %>% as.character()


PhysioData<-PhysioData %>% filter(BiopacSubject!="sp_9")
PhysioData<-PhysioData %>% filter(BiopacSubject!="sp_10")

#aggregate spielfeld 
PhysioData$Spielfeld<-as.numeric(
  as.character(
    unlist(PhysioData$Spielfeld)))%%16+1

#drop biopactrail (only from merging)

col <- c("Target","Measurement", "BiopacSubject","BiopacTrail",
         "Time","Gruppe","betweenCond","Spielfeld",
         "Factor1","Factor2")
PhysioData["Time"]<-as.numeric(as.character(
  unlist((PhysioData["Time"]))))
PhysioData[col]<-lapply(PhysioData[col],factor)
sapply(PhysioData, class)

if(F){  
tbo<-CreateTableOne(data = PhysioData )
summary(tbo) 
}
rm(col,PhysioData_org)
```

#### Plots Measurement by Subject und Ausschluss von VPn

```{r eval=F, include=T}
plotMeasurementBySubject<- function(measurement,string){
    measurement<-measurement %>% filter(Measurement==string)
    ggplot(data = drop_na(measurement),mapping = aes(x = BiopacSubject, 
                         y = Value, 
                         fill = Measurement)) +
                      geom_violin() +
                      geom_jitter(width = 0.2, alpha = 0.6) +
                      theme_classic()
}
#Phight:
PhysioData %>% filter(BiopacSubject!="nsp_14")  %>%  
  plotMeasurementBySubject("Phight")
#PRQ
PhysioData  %>% filter(BiopacSubject!="nsp_14") %>%  
  plotMeasurementBySubject("STev")
#QT
PhysioData %>% filter(BiopacSubject!="nsp_14") %>% 
  plotMeasurementBySubject("QT")
#QTwidth
PhysioData %>% filter(BiopacSubject!="nsp_14") %>%  
  plotMeasurementBySubject("QTwidth")
#RRi
PhysioData %>% filter(BiopacSubject!="nsp_14") %>%  
  plotMeasurementBySubject("RRi")
#Rhight
PhysioData %>% filter(BiopacSubject!="nsp_14") %>%  
  plotMeasurementBySubject("Rhight")
#ST
PhysioData %>% filter(BiopacSubject!="nsp_14") %>% 
  filter(Value<=3) %>%  plotMeasurementBySubject("ST")
#STev
PhysioData  %>% filter(BiopacSubject!="nsp_14") %>% 
  plotMeasurementBySubject("STev")

#exclude:

PhysioData_1 <- PhysioData %>% 
  filter(BiopacSubject!="nsp_14") %>% 
  filter(BiopacSubject!="nsp_3") %>% 
  filter(BiopacSubject!="nsp_7") %>% 
  filter(Measurement == "RRi") 
# Set 0 Values to NA -> 
# these trails were either not anwsered or the trail timedout -> 
#may lead to skrewed distribution!

PhysioData_1$VasSlide.RT<-if_else(
  PhysioData_1$VasSlide.RT!=0,PhysioData_1$VasSlide.RT,NULL)

saveRDS(PhysioData_1,"physio.RT")

```

Phight:
- nsp_14: grosse varianz nach oben - exclude
- nsp_17/20 bottem effekt - 0
- nsp_7: alles null

PRQ:
nicht interpretierbare Verteilungen

QT:
- nsp_14 riesige varianz exclude
- nsp_7 alles null

QTwidth:
sehr komische Verteilungen, viele nullen, bei 100 beschränkt ?? komisch

Rhight:
- nsp_3 viele 0 
- nsp_7 alles 0
- nsp_14 riesige varianz exclude
- nsp_4 grosse varianz

RRi:
- nsp_3 lot of values at 0 --> exclude
- nsp_14 grosse varianz

ST: 
- nsp_7 alles 0
- ausreisser bei nsp_16, nsp_4 und sp_14 (Werte zwischen 10 und 70) -> filtern
- nsp4,5,6 grosse varianz

STev:
- prakisch alle nsp haben bottem effekt
- nsp_4 grosse varianz
- nsp_14 riesen varianz exclude


--> nur analyse von RRi, da nur für dieses Mass eine Hypothese besteht. 
Aussschluss von VP: nsp_14,nsp_3,nsp_7

### TimePlot RRi and condition
```{r eval=F, include=T}
PhysioData_1<-read_rds("physio.RT")
 physio_RRi<-na.omit(filter(PhysioData_1,
                            (Measurement=="RRi")&(grepl("Stimulus",Target))
                            &(Time!=-2)))


qplot(Time, Value, data=physio_RRi, geom=c("boxplot"), 
   fill=Time, main="Kategorial Time (s) against BPM",
   xlab="Time(s)", ylab="BPM")


if (F){
  mg <- ggplot(physio_RRi, aes(x = VasSlide.RT, y =VasSlide.VAS ,
                               colour = factor(BiopacSubject)))
  + geom_point()
  mg + facet_grid(Factor1 + Factor2 ~ betweenCond)
}
```

### Behavior 
```{r eval=FALSE, include=T}
p_load(readxl)
sp<-read_excel(paste0(SysDir,
                      "SoOp/deprecated/Soccer deprecated/original/rawDaten/soccerPlayer/E-Prime Daten/Merge_VP01_VP30.xlsx")) 
sp$betweenCond<-"soccerPlayer"
nsp<-read_excel(paste0(SysDir,
                       "SoOp/deprecated/Soccer deprecated/original/rawDaten/nonSoccerPlayer/E-Prime Daten/Merge_VP01_VP30.xlsx")) 
nsp$betweenCond<-"nonSoccerPlayer"
behaviorData<-rbind(sp,nsp) %>% 
  select(c(Name,Gruppe,betweenCond,
           Spielfeld,VasSlide.RT,VasSlide.VAS))
behaviorData<-behaviorData %>% 
  filter(!grepl("Bsp",Spielfeld))
#create condition
behaviorData$Factor1<-if_else(
  grepl(paste(c("Selbst","Ingroup1","Ingroup2"),
              collapse="|"),behaviorData$Gruppe),"Us","Them")
behaviorData$Factor2<-if_else(
  grepl(paste(c("Selbst","Konkurrent"),collapse="|"),
        behaviorData$Gruppe),"Single",              
  if_else(grepl(paste(c("Ingroup1","Outgroup1"),
                      collapse="|"),behaviorData$Gruppe),
          "Team_1","Team_2")) 
#aggregate spielfeld 
behaviorData$Spielfeld<-as.numeric(
  as.character(unlist(behaviorData$Spielfeld)))%%16+1
behaviorData$VasSlide.RT<-if_else(
  behaviorData$VasSlide.RT!=0,behaviorData$VasSlide.RT,NULL)
#Modeling
saveRDS(behaviorData,"behavior.rds")
```

```{r include=FALSE, eval=F}
#plot
behaviorData<-read_rds("behavior.rds")
mg <- ggplot(behaviorData, 
             aes(x = VasSlide.RT, y =VasSlide.VAS ,
                 colour = factor(Name))) + geom_point()
mg + facet_grid(Factor1 + Factor2 ~ betweenCond)
```

# Baysian Modelling
Modelle müssen zuerst berechnet werden. Dazu bitte die entsprechenden R Scripte ausführen. 

#### Modeling Function (getBRMModel.R)
```{r eval=FALSE, code = readLines('mod/getBRMModel.R')}
```

#### Behavior (brms_b.R)
```{r eval=FALSE, code = readLines('mod/brms_b.R')}
```

#### Reaktion Times (brms_RT.R)

```{r eval=FALSE, code = readLines('mod/brms_RT.R')}
```

#### RRi (brms_RRi.R)

```{r eval=FALSE, code = readLines('mod/brms_RRi.R')}
```

### Behavior

```{r eval=F, include=T}
# EXECUTE brms_b.R to generate Models!

models<-data_frame(list("m0",
             "m1a",
             "m1b",
             
             "m2",
             "m2a",
             "m2b",
             "m2c"
             )) 
colnames(models)[1]<-"names"
models$measurement<-"VAS"
models <- split(models, seq(nrow(models)))
M<- lapply(models,
           function(x){read_rds(
             paste0("models/",x$measurement,
                    "_",x$names,".rds"))})

if (F){
  loo_behavior<- loo(M$`1`,M$`2`,M$`3`,M$`4`,M$`5`,M$`6`,M$`7`,
                     pointwise = F,
                     cores = parallel::detectCores ())
  saveRDS(loo_behavior,"models/loo_behavior.rds")
} else{
  loo_behavior<-read_rds("models/loo_behavior.rds")
}
lapply(M,function(x){x$formula})
loo_behavior
summary(M$`6`)
dat <- ggpredict(M$`6`,
                 terms = c("Factor1", "Factor2","betweenCond"),ppd=T)
saveRDS(dat,"Verhaltensdaten_plot.rds")
plot(dat,alpha = 0.05, dodge = 0.5)  + labs(title = "Verhaltensdaten")   +
  ylab("Visual Analoge Scale (1-100)") 
```

### Reaction Times
```{r eval=F, include=T}
# EXECUTE brms_RT.R to generate Models!
models<-data_frame(list("m0",
             "m1a",
             "m1b",
             
             "m2",
             "m2a",
             "m2b",
             "m2c"
             ))
colnames(models)[1]<-"names"
models$measurement<-"RT"
models <- split(models, seq(nrow(models)))
M_RT<- lapply(models,
              function(x){read_rds(
                paste0("models/",
                       x$measurement,"_",x$names,".rds"))})

if (F){
  loo_RT<- loo(M_RT$`1`,M_RT$`2`,M_RT$`3`,M_RT$`4`,
               M_RT$`5`,M_RT$`6`,M_RT$`7`,
               pointwise = F,
               cores = parallel::detectCores ())
  saveRDS(loo_RT,"models/loo_RT.rds")
} else{
  loo_RT<-read_rds("models/loo_RT.rds")
}
lapply(M_RT,function(x){x$formula})
loo_RT
summary(M_RT$`2`)

```


### RRi
```{r eval=F, include=T}
# EXECUTE brms_RRi.R to generate Models!
models<-data_frame(c(1:13))
colnames(models)[1]<-"names"
models$measurement<-"RRi"
models <- split(models, seq(nrow(models)))
M_RRi<- lapply(models,
               function(x){read_rds(
                 paste0("models/",
                        x$measurement,"_",
                        x$names,".rds"))})

if (F){
  loo_RRi<- loo(M_RRi$`1`,M_RRi$`2`,M_RRi$`3`,M_RRi$`4`,
                M_RRi$`5`,M_RRi$`6`,M_RRi$`7`,M_RRi$`8`,
                M_RRi$`9`,M_RRi$`10`,M_RRi$`11`,M_RRi$`12`,
                pointwise = F,cores = parallel::detectCores ()) 
  saveRDS(loo_RRi,"models/loo_RRi.rds")
} else {
  loo_RRi<-read_rds("models/loo_RRi.rds")
}
lapply(M_RRi,function(x){x$formula})
loo_RRi

print(M_RRi$`12`,digits=2)

p_load(ggeffects)

dat <- ggpredict(M_RRi$`12`, 
                 terms = c("Factor1", "Factor2","betweenCond"),ppd=T)
plot(dat,alpha = 0.05, dodge = 0.5 )+ labs(title = "Physiologie (RRi)")   +
  ylab("Herzrate in BPM") 
saveRDS(dat,"RRi_plot.rds")
read_rds("models/RRi_12.rds")
```