Load Libraries

library(ggplot2)
library(ggeffects)
library(psych)
library(tidyr)
library(lmerTest)
library(lme4)
library(foreign)
#library(readxl)

A. Gather Psychophys Data

#subdirs <- list.dirs(path = "/Volumes/Groups/AHeller_Lab/Undergrad/HealthyUStudy/Session 7/Psychphysiology/Filtered_converted", recursive = F)

B. Pull specific subject data

#summarydat<-matrix(0,nrow = (length(subdirs)), ncol = 4)

# for (i in 1:39){    #i =24 & 27 is a manual one (5320) ||| i=27,
#   
#   thissubdir <- subdirs[i]
#   thissub <- substr(thissubdir, 98,101)
#   thissubdir <- paste(thissubdir, "/output data",sep="")
#   #grab just the hrv output files
#   subfiles <- list.files(path = thissubdir, pattern = "HRV", recursive = T)
#       
#   if (thissub == "5300"){#for 5300,5301,5304
#       subfiles <- list.files(path = thissubdir, pattern = "PM", recursive = T)
#   }
#   if (thissub == "5301"){#for 5300,5301,5304
#       subfiles <- list.files(path = thissubdir, pattern = "PM", recursive = T)
#   }
#   if (thissub == "5304"){#for 5300,5301,5304
#       subfiles <- list.files(path = thissubdir, pattern = "PM", recursive = T)
#   }
#   if (thissub == "5305"){#
#       subfiles <- list.files(path = thissubdir, recursive = T)
#       subfiles <- subfiles [c(26:30)]
#   }
#    if (thissub == "5306"){
#       subfiles <- list.files(path = thissubdir, pattern = "hrv", recursive = T)
#   }
#   if (thissub == "5203"){#
#       subfiles <- subfiles[c(1:5)]
#   }
#   if (thissub == "5320"){#
#       subfiles <- list.files(path = thissubdir, recursive = T)
#   }
#   if (thissub == "5320"){#
#       subfiles <- list.files(path = thissubdir, recursive = T)
#   }
#   
#   #in case they are saved lower case...
#   if (length(subfiles)==0) {
#       subfiles <- list.files(path = thissubdir, pattern = "hrv", recursive = T)
#   }
# 
#       #see/save what we are working with
#       print(paste(".....", thissub, "has a total of", length(subfiles), "runs to read in ....."))
#       summarydat[i,1] <- thissub
#       summarydat[i,2] <- length(subfiles)
#       
#       for (r in 1:3){
#     
#       runtext <- paste("run",r,sep = "")
#       filenum <- grep(runtext, subfiles)
#       #filenum <- r
#       #hopefully this with stop it from breaking...
#       if (length(filenum)!=0) {
#         filename <- paste(thissubdir,"/",subfiles[filenum],sep="")
#         filedat <- read_xlsx(filename, col_names = F)
#         filedat_t <- as.data.frame(t(filedat))
# 
#         red_dat <- as.data.frame(filedat_t[,c(1,3,5,6,9:12)])
#         red_dat <- red_dat[-c(1),]
#         colnames(red_dat) <- c("SegmentNumber", "StartTime", "EndTime", "SegmentDuration", "MeanHeartRate", "RSA", "MeanIBI", "NumberRsFound")
# 
#         red_dat$id <- thissub
#         red_dat$run <- r
#         
#           #if (r==1){
#           #summarydat[i,3] <- length(red_dat$SegmentNumber)
#           #summarydat[i,4] <- as.character(red_dat$SegmentDuration[1])
#           #}
#         
#         red_dat <- red_dat[c(9,10,1:8)]
#         red_dat <- red_dat[-c(12),]
#         #alldatmat <- red_dat #just for first sub, first run
#         alldatmat <- rbind(alldatmat,red_dat)
#         } #close if/filenum
# 
#     } #close run loop
#   } #close subject loop
# 
# #write.csv(alldatmat, "/Volumes/Groups/AHeller_Lab/Undergrad/HealthyUStudy/Session 7/Psychphysiology/Summary_051022.csv")
# #write.csv(alldatmat, "/Users/nikki/Desktop/Research/Dissertation/Analyses/Physio/Summary_051022.csv")
# write.csv(alldatmat, "/Users/nikki/Desktop/Research/Dissertation/Analyses/Physio/FullSample_Summary_071522.csv")

alldatmat <- read.csv("/Users/nikki/Desktop/Research/Dissertation/Analyses/Physio/FullSample_Summary_071522.csv")
alldatmat <- alldatmat[-c(1)]
alldatmat$SegmentNumber <- as.numeric(alldatmat$SegmentNumber)
alldatmat_segs <- subset(alldatmat,alldatmat$SegmentNumber < 12)

C. Format trial information to line up with phys data

scandat <- read.csv("/Users/nikki/Desktop/Research/HealthyU_Scanning_Local/Emily_HUresources/EventTimingData_long_EM.csv")

scandat$Trial <- NA
scandat$Trial <- rowSums(scandat[c(9:13)], na.rm = T) #create a single variable with the run-wise trial numbers
scandat <- scandat[c(1,7,26,14:20,22,23)] #reduct to need vars 
scandat$Block <- as.factor(scandat$Block)
scandat$Block <- as.numeric(scandat$Block)

#finally, reshape so that the trial types (event, feel and int are ~wide~) I have already done this somewhere
  # The arguments to spread():
  # - data: Data object
  # - key: Name of column containing the new column names
  # - value: Name of column containing values

scandat <- scandat[-c(7,8,10:12)]
scandat_wide <- spread(scandat, DataType, Response)
scandat_wide <- scandat_wide[-c(6)]
scandat_wide$EventType <- as.character(scandat_wide$EventType)
subjects <- unique(alldatmat$id)
# alldatmat$EventType <- NA
# alldatmat$FeelingRate <- NA
# alldatmat$IntensityRate <- NA 
# alldatmat$subcheck <- NA
# alldatmat$blockcheck <- NA
# alldatmat$trialcheck <- NA
runs <- c(1:5)

colnames(scandat_wide)
## [1] "Subject"       "Block"         "Trial"         "Event"        
## [5] "EventType"     "FeelingRate"   "IntensityRate"
colnames(scandat_wide) <- c("id","run","SegmentNumber","Event","EventType","FeelingRate","IntensityRate")
mergetest <- merge(alldatmat, scandat_wide, by = c("id","run","SegmentNumber"))

# for (s in 1:length(subjects)){
#   subjecttrials <- which(scandat_wide$Subject == subjects[s]) 
#   subjectdata <- scandat_wide[subjecttrials,]
#   match_trials <- which(alldatmat$id == subjects[s])
#       for (r in 1:length(runs)){
#       runtrials <- which(subjectdata$Block == runs[r])
#       match_runtrials <- which(alldatmat$run[match_trials] == runs[r])
#       match_runtrials <- match_trials[match_runtrials]
#           if (length(match_runtrials) == 11){
#             alldatmat[match_runtrials,c(11:16)] <- subjectdata[runtrials,c(5:7,1:3)]
#           } else if (length(match_runtrials) == 10) {
#             newruntrials <- c(runtrials[2]:runtrials[11])
#             alldatmat[match_runtrials,c(11:16)] <- subjectdata[newruntrials,c(5:7,1:3)]
#           } else if (length(match_runtrials) == 12) {
#             newmatch_runtrials <- c(match_runtrials[1]:match_runtrials[11])
#             alldatmat[newmatch_runtrials,c(11:16)] <- subjectdata[runtrials,c(5:7,1:3)]
#           } else 
#             alldatmat[match_runtrials,c(11:16)] <- NA
#        }
# }

#write.csv(mergetest,"/Volumes/Groups/AHeller_Lab/Undergrad/HealthyUStudy/Session 7/Psychphysiology/Summary_wEventLabels_053122.csv" )
#write.csv(mergetest,"/Users/nikki/Desktop/Research/Dissertation/Analyses/Physio/Summary_wEventLabels_060822_5329.csv")
write.csv(mergetest,"/Users/nikki/Desktop/Research/Dissertation/Analyses/Physio/Summary_wEventLabels_071522_FULL.csv")

alldatmat_b <- alldatmat
alldatmat <- mergetest

alldatmat$EventType <- as.factor(alldatmat$EventType)

alldatmat <- alldatmat[ order( alldatmat$id, alldatmat$run, alldatmat$SegmentNumber),]

D. Multilevel model to test whether HR changes from condition to condition

OR by emotional intensity

Formatting

### PREP 
#alldatmat<-alldatmat_clean
alldatmat$id<-as.numeric(alldatmat$id)
alldatmat$RSA<-as.numeric(alldatmat$RSA)
alldatmat$MeanHeartRate<-as.numeric(alldatmat$MeanHeartRate)

RSAgrpmeans<-aggregate(alldatmat$RSA,list(alldatmat$id),mean,na.rm=TRUE)
names(RSAgrpmeans)<-c("id","prsmn.RSA")
alldatmat<-merge(alldatmat,RSAgrpmeans,by="id") #merge group means into original data
alldatmat$cwp.RSA<-alldatmat$RSA-alldatmat$prsmn.RSA # person centered

HRgrpmeans<-aggregate(alldatmat$MeanHeartRate,list(alldatmat$id),mean,na.rm=TRUE)
names(HRgrpmeans)<-c("id","prsmn.HR")
alldatmat<-merge(alldatmat,HRgrpmeans,by="id") #merge group means into original data
alldatmat$cwp.HR<-alldatmat$MeanHeartRate-alldatmat$prsmn.HR # person centered

Intgrpmeans<-aggregate(alldatmat$IntensityRate,list(alldatmat$id),mean,na.rm=TRUE)
names(Intgrpmeans)<-c("id","prsmn.Int")
alldatmat<-merge(alldatmat,Intgrpmeans,by="id") #merge group means into original data
alldatmat$cwp.Int<-alldatmat$IntensityRate-alldatmat$prsmn.Int # person centered

Modeling

### MODELS
  
  # 1. HR ~ Event is not sig
  HR_event_model <- lmer(MeanHeartRate ~ EventType + (1|id), REML = FALSE, data = alldatmat)
  summary(HR_event_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: MeanHeartRate ~ EventType + (1 | id)
##    Data: alldatmat
## 
##      AIC      BIC   logLik deviance df.resid 
##  11252.0  11280.1  -5621.0  11242.0     2040 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.7885 -0.5495 -0.0861  0.4429  9.3289 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 85.82    9.264   
##  Residual             12.81    3.579   
## Number of obs: 2045, groups:  id, 38
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      67.1738     1.5105   38.5680  44.472   <2e-16 ***
## EventTyperum      0.1282     0.2005 2007.0047   0.639    0.523    
## EventTypeworry    0.2705     0.2005 2007.0047   1.349    0.177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) EvntTypr
## EventTyperm -0.076         
## EvntTypwrry -0.076  0.572
  # 2. HR ~ Int is sig!
  HR_intensity_model <- lmer(MeanHeartRate ~ IntensityRate + (1 |id), REML = FALSE, data = alldatmat)
  summary(HR_intensity_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: MeanHeartRate ~ IntensityRate + (1 | id)
##    Data: alldatmat
## 
##      AIC      BIC   logLik deviance df.resid 
##  11117.8  11140.3  -5554.9  11109.8     2017 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.7367 -0.5501 -0.0863  0.4525  9.2647 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 85.82    9.264   
##  Residual             12.79    3.577   
## Number of obs: 2021, groups:  id, 38
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.677e+01  1.516e+00 3.914e+01  44.046  < 2e-16 ***
## IntensityRate 2.621e-01  8.475e-02 1.985e+03   3.093  0.00201 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## IntensityRt -0.121
  # 3. Decomposing win/bw
  HR_intensity_model_Dec <- lmer(MeanHeartRate ~ cwp.Int + prsmn.Int + (1 |id), REML = FALSE, data = alldatmat)
  summary(HR_intensity_model_Dec)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: MeanHeartRate ~ cwp.Int + prsmn.Int + (1 | id)
##    Data: alldatmat
## 
##      AIC      BIC   logLik deviance df.resid 
##  11119.7  11147.7  -5554.8  11109.7     2016 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.7369 -0.5502 -0.0866  0.4527  9.2646 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 85.55    9.250   
##  Residual             12.79    3.577   
## Number of obs: 2021, groups:  id, 38
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   69.81034    9.00342   38.00013   7.754 2.41e-09 ***
## cwp.Int        0.26273    0.08477 1983.00461   3.099  0.00197 ** 
## prsmn.Int     -1.14400    4.11145   37.99905  -0.278  0.78233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) cwp.In
## cwp.Int    0.000       
## prsmn.Int -0.986  0.000

Plotting

### PLOT

myplot <- ggplot(data = alldatmat, aes(x= cwp.Int, y = cwp.HR, group = id, color = id)) +
        geom_line(stat="smooth",method = "lm", size = .5, alpha = 0.5) +  
        geom_smooth(aes(group=1), color="black", alpha = .4, size = 1.5, method = "lm", se = F) +
        #scale_color_gradient(low="coral2", high="purple") + 
        xlab("Within-person Scanner Intensity Ratings") +
        ylab("Within-person Trial-wise Heart Rate") 
myplot + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))