Commit 5a58744d authored by Amit Goldenberg's avatar Amit Goldenberg

a few updates to the code

parent b34e7a8b
......@@ -19,9 +19,15 @@ We assume that these features of individual-social mismatch trigger two differen
```{r PREPARE R, include=FALSE}
if (!require("pacman")) install.packages("pacman")
pacman::p_load("ggplot2","lme4",'lmerTest','dplyr',"tidyr",'Rmisc','GGally',"RColorBrewer", "sjPlot", "MuMIn")
# if (!require("pacman")) install.packages("pacman")
# pacman::p_load("ggplot2","lme4",'lmerTest','dplyr',"tidyr",'Rmisc','GGally',"RColorBrewer", "sjPlot", "MuMIn")
*Amit*- I had issues with pacman
library (ggplot2); library(lme4) ; library (lmerTest) ; library (dplyr); library (tidyr); library(Rmisc);
library (GGally); library (RColorBrewer); library (sjPlot)
# make plotting theme
pl.pd <- position_dodge(0.2); pl.sc <- 2.5; pl.trsp = 0.6;
pal <- brewer.pal("Blues", n=9); pl.bkg <- pal[2]; pl.mid <- pal[5]; pl.foc <- pal[8]
......@@ -44,53 +50,54 @@ We assume that these features of individual-social mismatch trigger two differen
legend.title = element_blank(),legend.position = "right"
)
```
```{r IMPORT DATA, include=FALSE}
#Amit's folder
# setwd("C:/Users/Amit/Dropbox/Research/emotional conformity-non conformity/motivated_contation_EEG_2017/Results")
setwd("C:/Users/Amit/Dropbox/Research/emotional conformity-non conformity/motivated_contation_EEG_2017/Results")
#andero's folder
setwd("C:/Users/ro/OneDrive - Tartu Ülikool/Teadus/REGULATION/16 motivated contagion/shared/")
#setwd("C:/Users/ro/OneDrive - Tartu Ülikool/Teadus/REGULATION/16 motivated contagion/shared/")
## 1. PERPARE DATA
# import single trial behavioral + EEG file where each row = trial
d <- read.csv2("eeg_data_motivated_contagion_revised version2.9.17.csv")
head(d)
d <- read.csv2("eeg_data_motivated_contagion_revised version2.9.17.csv")
head(d)
# remove group ratings from "none" conditions
d$grate = ifelse(d$grate == 0, NA, d$grate)
# remove group ratings from "none" conditions
d$grate = ifelse(d$grate == 0, NA, d$grate)
# remove strange trials
d <- d[d$rate %in% c(1:7),]
d$rt <- as.numeric(as.character(d$rt))
head(d)
# remove strange trials
d <- d[d$rate %in% c(1:7),]
d$rt <- as.numeric(as.character(d$rt))
head(d)
# remove the electrode location labels from the colnames
colnames(d) <- as.character(lapply(colnames(d), function(x) strsplit(x,"_")[[1]][1]))
## 2. MAKE ROW = MEASUREMENT FILE
# make the files longer (we have several time windoes from each trial in separate columns, we want create a file where each row has a single DV
# remove the electrode location labels from the colnames
colnames(d) <- as.character(lapply(colnames(d), function(x) strsplit(x,"_")[[1]][1]))
head(d)
# 2. MAKE ROW = MEASUREMENT FILE
# make the files longer (we have several time windoes from each trial in separate columns, we want create a file where each row has a single DV
# base is the information about what's going on each trial
base <- d[,c(1:which(colnames(d)=="Subject"))]
# we do the reshapting separately for each type of EEG correlate - ERP or specific frequency bands.
alpha <- merge(base,gather(d,window, alpha, alpha1:alpha4)[,c("Subject", "urevent", "window", "alpha")],c("Subject", "urevent"))
theta <- merge(base,gather(d,window, theta, theta1:theta4)[,c("Subject", "urevent", "window", "theta")],c("Subject", "urevent"))
delta <- merge(base,gather(d,window, delta, delta1:delta4)[,c("Subject", "urevent", "window", "delta")],c("Subject", "urevent"))
alpha$window <- gsub("alpha", "", alpha$window); theta$window <- gsub("theta", "", theta$window); delta$window <- gsub("delta", "", delta$window)
# base is the information about what's going on each trial
base <- d[,c(1:which(colnames(d)=="Subject"))]
# we do the reshapting separately for each type of EEG correlate - ERP or specific frequency bands.
alpha <- merge(base,gather(d,window, alpha, alpha1:alpha4)[,c("Subject", "urevent", "window", "alpha")],c("Subject", "urevent"))
theta <- merge(base,gather(d,window, theta, theta1:theta4)[,c("Subject", "urevent", "window", "theta")],c("Subject", "urevent"))
delta <- merge(base,gather(d,window, delta, delta1:delta4)[,c("Subject", "urevent", "window", "delta")],c("Subject", "urevent"))
alpha$window <- gsub("alpha", "", alpha$window); theta$window <- gsub("theta", "", theta$window); delta$window <- gsub("delta", "", delta$window)
# since all frequency bands have the same number of time windoes, we concatenate all frequency data into a single frame
ersp <- merge(merge(alpha,theta[,c("Subject", "urevent", "window", "theta")],c("Subject", "urevent", "window")),
delta[,c("Subject", "urevent", "window", "delta")],c("Subject", "urevent", "window"))
erp <- merge(base,gather(d,window, erp, P2:LPP4)[,c("Subject", "urevent", "window", "erp")],c("Subject", "urevent"))
ersp <- merge(merge(alpha,theta[,c("Subject", "urevent", "window", "theta")],c("Subject", "urevent", "window")),delta[,c("Subject", "urevent", "window", "delta")],c("Subject", "urevent", "window"))
erp <- merge(base,gather(d,window, erp, P2:LPP4)[,c("Subject", "urevent", "window", "erp")],c("Subject", "urevent"))
# add stuff to base just so the following loop would work
base$window <- 1
base$value <- 0
rm(alpha,theta,delta)
# separate phases (we now make the file shorter by bringing the individual and social phase into the same row)
for(file in c("base", "ersp", "erp")){
# separate phases (we now make the file shorter by bringing the individual and social phase into the same row)
for(file in c("base", "ersp", "erp")){
eval(parse(text=paste0("d <- ",file)))
......@@ -167,6 +174,10 @@ head(erp)
colnames(d)
# we have to remove situations in which the first rating is 8
d = subset (d, rate.indiv < 8 & rate.indiv >1 )
erp = subset (erp, rate.indiv < 8 & rate.indiv >1 )
base = subset (base, rate.indiv < 8 & rate.indiv >1 )
```
......@@ -259,7 +270,7 @@ ggplot(base,aes(x=rate.indiv, y = grate, color = condition)) + geom_jitter() + g
We kind of have a dilemma. On the one hand, we could stick to analyzing the conditions. On the other hand, each conditions contains quite a bit of variance so we could also treat the group rating as a continuous variable. The added benefit of the continuous approach is that we don't have to worry about if the conditions are perfectly balanced after we've lost trials due to EEG preprocessing.
#amit's comment: I vote for grate as a countious.
*amit*: I vote for grate as a countious.
Now, if we wanted to use continuous independent variables, which ones make the most sense?
......@@ -275,6 +286,15 @@ ggpairs(base[,c("grate", "rate.social", "rate.indiv")]) + pl.theme + labs(title
I did this by predicting second phase ratings for first phase ratings in the none conditions using a GLMM with variable intercepts and slopes. I then took the intercepts and sloped for each subject and computed expected second phase ratings for each picture.
*Amit* - I like this direction idea but it has a few isssues.
1. We don't really know what is the psychologocial process that participants go through during this stage. Second
2. The none conditions just come up once in 4 trials. This means that the habituation coef is not a per-trial coef but per 4 tirals.
3. There are a few types of habituation in this file: (1)between the first phase and the second phase. (2) whitin each phase. You are suggesting to deal just with one.
Here are a few alternatives:
1. I vote to start with just a difference score and contorl for trial number.
2. If we want to take an individual level habituation coef at each phase we should just look at change per trial for all conditions.
Here're the predicted ratings together with actual ones.
```{r, warning=F}
ggpairs(base[,c("rate.indiv", "rate.social", "rate.predict")]) + pl.theme + labs(title = "Observed and predicted invidual ratings") + pl.theme
......@@ -290,6 +310,10 @@ ggplot(base,aes(x=mismatch2, y = rate.social, color = condition)) + geom_jitter(
I operationalized the direction for now just based on quantiles of the mismatch variable - below 1/3 = -1; above 2/3 = 1; and the rest = 0.
*amit* - remember we have to remove the first phase situations in which the first rating is either 1 or 8.
```{r, warning=F}
ggpairs(base[,c("mismatch2", "direction2", "magnitude2")]) + pl.theme + labs(title = "Analytic variables") + pl.theme
```
......@@ -299,19 +323,33 @@ We can use direction and magnitude as the main independent variables and go and
## Behavior
First, let's apply this approach to behavioral data (ratings from the second phase excluding none trials).
```{r, warning=F, results="hide", message=F}
head(base)
base$difNum = base$rate.social - base$rate.indiv
base$conCont = base$grate- base$rate.indiv
fit <- lmer(rate.social ~ rate.predict + direction2 + magnitude2 + (1|order.social) + (1|Subject), base[base$condition!="none",])
summary(fit)
anova(fit)
sjp.lmer(fit, type = "fe.std", y.offset = .3, prnt.plot = F)$plot + pl.theme +
labs(subtitle = paste0("Model explains ",round(r.squaredGLMM(fit)[[2]]*100,1),
"%, fixed effects explain ", round(r.squaredGLMM(fit)[[1]]*100,1),"%."),
title = "Standradized effect sizes for main model")
sjp.lmer(fit, type = "eff", show.ci = T, prnt.plot = F)$plot + pl.theme + ggtitle("Marginal predictor effects")
### amit
fit = lmer (scale(difNum, scale =F) ~ condition+order.indiv -1 + (1|photo)+ (1|Subject), base[base$condition!="none",]) ; summary (fit)
fit = lmer (rate.social ~ condition +rate.indiv+ (1|photo)+ (1|Subject), base[base$condition!="none",]) ; summary (fit)
head(base)
```
......@@ -329,7 +367,16 @@ ggplot(base,aes(x=mismatch, y = rate.dif)) + stat_summary(fun.data = "mean_cl_bo
Given that we are interested in the time course of effects, I've exported consequitive time-windows (P2,P3, and 4 LPPs with 250 ms increments). The time windows is entered as a factor to the GLMM and we expect it to show interactions whenever time matters.
```{r, warning=F, results="hide", message=F}
head (erp)
erp$difNum = erp$rate.social - erp$rate.indiv
erp$conCont = erp$grate- erp$rate.indiv
erp$erp.dif = erp$erp.social - erp$erp.indiv
fit <- lmer(erp.social ~ -1 + rate.predict*window + direction2*window + magnitude2*window + (1|order.social) + (1|Subject), erp[erp$condition!="none",])
#if you are removing the intecept you have to center the outcome variable.
anova(fit)
sjp.lmer(fit, type = "fe.std", y.offset = .3, prnt.plot = F)$plot + pl.theme +
labs(subtitle = paste0("Model explains ",round(r.squaredGLMM(fit)[[2]]*100,1),
......@@ -341,6 +388,97 @@ ggplot(erp,aes(x=mismatch, y = erp.dif)) + stat_summary(fun.data = "mean_cl_boot
##amit - analyze per window
#let's first see in which windows erp match each other
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.social ~ erp.indiv +(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5] ,group = .[[1]] )) ;res# th
#p1 - yes, p2 - yes, lpp1 - yes, lpp2 - yes, lpp3 - no , lpp4 - no
##amit - analyze per window
#let's first see in which windows erp match each other
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.social ~ erp.indiv+rate.indiv+rate.social +(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5] ,group = .[[1]] )) ;res# th
#let's see in which windews erp match ratings
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.social ~ rate.social +(1|photo)+(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5] ,group = .[[1]] )) ;res# th
# this is significant in all windows but the peak is lpp1 - Let's look at the individual rating vs lpp
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(rate.social ~ rate.indiv +(1|photo)+(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5] ,group = .[[1]] )) ;res# th
# this is significant in all windows but the peak is lpp1 - Let's look at the individual rating vs lpp
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.indiv ~ rate.indiv +(1|photo)+(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5] ,group = .[[1]] )) ;res# th
#interstingly the peak here comes later - at lpp2
#let's see in which windews erp match ratings
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.social ~ grate+ rate.social+rate.indiv+erp.indiv +(1|Subject), data = .))))
res=r %>%
do(data.frame (coef_grate = summary(.$r)$coefficients[2], p_grate =summary(.$r)$coefficients[2,5],coef_social = summary(.$r)$coefficients[3],p_social =summary(.$r)$coefficients[3,5] ,group = .[[1]] )) ;res# th
# it seems that the peak is again in lpp1
## let's create a difference score
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.social ~ grate+ rate.social+rate.indiv+erp.indiv +(1|photo)+(1|Subject), data = .))))
res=r %>%
do(data.frame (coef_grate = summary(.$r)$coefficients[2], p_grate =summary(.$r)$coefficients[2,5],coef_social = summary(.$r)$coefficients[3],p_social =summary(.$r)$coefficients[3,5] ,group = .[[1]] )) ;res# th
# it seems that the peak is again in lpp1
head(erp)
######## let's look at difference in ERP###
r= erp[erp$condition!="none",] %>%
group_by(window) %>%
do(r = print(summary(lmer(erp.dif ~ conCont +(1|photo)+(1|Subject), data = .))))
res=r %>%
do(data.frame (coef = summary(.$r)$coefficients[2], p =summary(.$r)$coefficients[2,5],group = .[[1]] )) ;res# th
# it seems that the peak is again in lpp1
head(erp)
```
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment