Commit 423264a6 authored by Amit Goldenberg's avatar Amit Goldenberg

added files from Daniel's computer

parent 095d54bd
%% EEG ANALYSES FOR AMIT'S MOTIVATED CONTAGION EXPERIMENT %%
clear all; close all; clc;
%% SETTINGS %%
do = struct('import',0, 'reimport',0, ... % flow control - a structure that monitors our progress
'ica',0, 'reica',0, ...
'epoch',1, 'reepoch',0, ...
do = struct('import',0 , 'reimport',0, ... % flow control - a structure that monitors our progress
'ica',0, 'reica',0, ... % rieca is an instruction whether to override ICA same with the others
'epoch',0, 'reepoch',0, ...
'build', 1, 'analyze', 0);
% Andero paths (change)
%dataDir = 'L:\REGULATION\Amit 2016\'; % directory for large datafiles
%scriptDir = 'S:\Teadus\affect regulation\16 motivated contagion\'; % directory for scripts and other files (mirrored to gitLab)
%eeglabDir = 'C:\Program Files\EegLab\eeglab\'; % eeglab program directory
% dataDir = 'L:\REGULATION\Amit 2016\'; % directory for large datafiles
% scriptDir = 'S:\Teadus\affect regulation\16 motivated contagion\'; % directory for scripts and other files (mirrored to gitLab)
% eeglabDir = 'C:\Program Files\EegLab\eeglab\'; % eeglab program directory
%Amit paths
% dataDir = 'C:\Users\Amit\OneDrive\Study_11_EEG_pilot\test\' % directory for large datafiles
% scriptDir = 'C:\Users\Amit\OneDrive\Study_11_EEG_pilot\Results\' % directory for scripts and other files (mirrored to gitLab)
% eeglabDir = 'C:\Program Files\MATLAB\R2016a\toolbox\eeglab13_5_4b\'; % eeglab program directory. I am not sure I need this if I comment
% Amit paths
% dataDir = 'C:\Users\Amit\OneDrive\Study_11_EEG_pilot\test\'
% scriptDir = 'C:\Users\Amit\OneDrive\Study_11_EEG_pilot\Results\'
% eeglabDir = 'C:\Program Files\MATLAB\R2016a\toolbox\eeglab13_5_4b\';
%Daniel paths
dataDir = 'D:\emotion contagion eeg\amit_conformity_data' % directory for large datafiles
scriptDir = 'D:\emotion contagion eeg\amit_conformity\EEG_analysis_conformity_task' % directory for scripts and other files (mirrored to gitLab)
eeglabDir = 'C:\Program Files\MATLAB\R2016a\toolbox\eeglab_current\eeglab13_4_4b\'; % eeglab program directory. I am not sure I need this if I comment
% Daniel paths
dataDir = 'D:\emotion contagion eeg\amit_conformity_data\';
scriptDir = 'D:\emotion contagion eeg\amit_conformity\EEG_analysis_conformity_task\';
eeglabDir = 'C:\Program Files\MATLAB\R2016a\toolbox\eeglab_current\eeglab13_4_4b\';
% RELATIVE PATHS (don't change)
addpath(scriptDir); %addpath(eeglabDir); % add the second part if it doesn't wrok % add scripts and eeglab to MATLAB path for this session
impdir = [dataDir, 'raw\']; % raw EEG files
indir = [dataDir,'import\']; mkdir(indir); % EEG files converted to EEGLAB format
icadir = [dataDir,'output\ica\']; mkdir(icadir); % EEG files with ICA solutions
epochdir = [dataDir,'output\epoch\']; mkdir(epochdir); % fully preprocessed EEG files
addpath(scriptDir); %addpath(eeglabDir); % add the second part if it doesn't wrok % add scripts and eeglab to MATLAB path for this session
impdir = [dataDir, 'raw\']; % raw EEG files
indir = [dataDir,'import\']; mkdir(indir); % EEG files converted to EEGLAB format
icadir = [dataDir,'output\ica\']; mkdir(icadir); % EEG files with ICA solutions
epochdir = [dataDir,'output\epoch_f01\']; mkdir(epochdir); % fully preprocessed EEG files
% DESIGN
d.photoType = {'emo' 'ntr'};
d.phase = {'indiv' 'social'};
d.condition = {'low' 'high' 'same' 'none'};
eeglab; % open eeglab
%% IMPORT %%
......@@ -47,6 +53,8 @@ while do.import==1; % fl
% indexing trialinfo with a subset of trialinfo. 2:end - rows I want to
% take. strcom - look for string call screen in the first row.
subList = unique([trialinfo{2:end,strcmp('student_id', trialinfo(1,:))}]); % build a cell with unique subject IDs
subList = setdiff(subList,[6065275 6002886 5934980]); % leave out subject with trigger problems
% unique just filters our repeats
for s = subList; % loop through subject IDs
if ~do.reimport; if any(strcmp([num2str(s),'.set'], {inlist(:).name})); continue; end; end % if a given subject has already been imported, skip if required in do
......@@ -71,19 +79,14 @@ while do.import==1; % fl
EEG = pop_resample( EEG, 500); % downsample to 500 Hz
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
% RE-REFERENCE
EEG = pop_reref(EEG, [23 24] ,'refloc',...
struct('labels',{'Cz'},'Y',{0},'X',{5.4492e-16},'Z',{8.8992},'sph_theta',{0},'sph_phi',{90},'sph_radius',{8.8992},'theta',{0},'radius',{0},'type',{'REF'},'ref',{''},'urchan',{[]},'datachan',{0}));
% DEFINE CHANNEL LOCATIONS
labels = {'Fp1' 'Fp2' 'F3' 'F4' 'C3' 'C4' 'P3' 'P4' 'O1' 'O2' 'F7' 'F8' 'T3' 'T4' 'T5' 'T6' 'Fz' 'Cz' 'Pz' 'Oz' ...
'21' '22' 'LM' 'RM' '25' '26' '27' 'FCz' 'EOG1' 'EOG2' 'EOG3' 'EOG4'};
% create a vector of channel labels in the order of channels in the data
for l = 1:length(labels); EEG.chanlocs(l).labels = labels{l}; end % write the new labels to the EEGLAB file
% IMPORT CHANNEL LOCATIONS
EEG=pop_chanedit(EEG, 'load',{[scriptDir 'locfile.ced'] 'filetype' 'autodetect'});
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
% RE-REFERENCE
[a, refChans] = intersect({EEG.chanlocs.labels},{'LM' 'RM'}); % find and reference to average mastoids
EEG = pop_reref(EEG, refChans, 'keepref', 'off');
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
% ADD EVENT INFORMATION
EEG.subject = num2str(s); % subject name
subinfo = trialinfo(1+find([trialinfo{2:end,strcmp('student_id',trialinfo(1,:))}]==s),:); % make a headerless matrix of trialinfo from the selected subject
......@@ -109,15 +112,22 @@ while do.import==1; % fl
row = find([subinfo{:,strcmp('order2',trialinfo(1,:))}]==count.two);
if isempty(row); continue; end % if the trial number is missing from behavioural file, continue
row = row(1); % if there are more than one trial numbers in behavioural file, take the first
for f = {'student_id', 'date', 'id', 'condition', 'photo', 'photoType'};
for f = {'student_id', 'date', 'id', 'condition', 'photo', 'photoType', 'grate'};
[EEG.event(e:e+2).(f{:})] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))});
end
[EEG.event(e:e+2).phase] = deal(2);
for f = {'rate2','rt2', 'picot2', 'rateot2', 'order2', 'groupot2'};
[EEG.event(e:e+2).(f{:}(1:end-1))] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))});
end
end
end
end
if any(strcmp({'PIC1' 'PIC3'}, EEG.event(e).type))
[EEG.event(e:e+2).photoType] = deal(d.photoType{EEG.event(e).photoType});
[EEG.event(e:e+2).phase] = deal(d.phase{EEG.event(e).phase});
if any(ismember(EEG.event(e).condition,[1:4]))
[EEG.event(e:e+2).condition] = deal(d.condition{EEG.event(e).condition});
end
end
end
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); % save changes to EEGLAB structure within Matlab
EEG = pop_saveset(EEG, 'filename', EEG.subject, 'filepath', indir, 'savemode', 'twofiles'); % save changes to the import directory on disk
end
......@@ -130,6 +140,7 @@ while do.ica == 1
% class of algo that seperates the variances into max independent sources
inlist = dir([indir , '*.set']); icalist = dir([icadir, '*.set']); %create a list of all the .set files in thie folder % update file lists
for s = 1:length(inlist); %we are looping through the list
if ~do.reica; if any(strmatch(inlist(s).name, {icalist(:).name})); continue; end; end
% if don't want to overwrite (this is a just a loop to get you
% started from the middle in cases your computer crashed) and if
......@@ -150,7 +161,7 @@ while do.ica == 1
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
[EEG ALLEEG(1).rejhist.preICA] = pop_rejspec(EEG, 1, [1:EEG.nbchan], -35, 35, 15, 30, 0, 0); %cleaning the data of muscle noise % remove segments with high muscle noise (15-30 Hz)
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);% whenever we see a lot of power for specific point.
EEG = pop_runica(EEG, 'extended', 1, 'interupt', 'on'); %actually running the ICA % find an ICA solution
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
......@@ -169,14 +180,30 @@ while do.epoch == 1;
icalist = dir([icadir,'*.set']); epochlist = dir([epochdir,'*.set']); % refresh different folder content lists
inlist = dir([indir,'*.set']);
for s = 1:length(icalist); % we are looking through the ica directory
if ~do.reepoch; if any(strmatch(icalist(s).name, {epochlist(:).name})); continue; end; end % overwrite check
ALLEEG = []; EEG = []; CURRENTSET = 1; % clean workspace and load dataset from ICA folder
% temporary event fix
EEG = pop_loadset('filename', icalist(s).name, 'filepath', indir);
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
EEG = pop_loadset('filename', icalist(s).name, 'filepath', icadir);
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET+1);
EEG.event = ALLEEG(1).event;
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
EEG = pop_saveset(EEG, 'filename', EEG.filename, 'filepath', icadir, 'savemode', 'twofiles'); % save the update ICA file back to disk
ALLEEG = []; EEG = []; % clean workspace and load dataset from ICA folder
EEG = pop_loadset('filename', icalist(s).name, 'filepath', icadir);
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
if isempty(find(EEG.reject.gcompreject)); % if the file does not have markings for ocular ICs ...
EEG = pop_epoch(EEG, {'CROS'}, [0 4], 'epochinfo', 'yes'); % epoch just for visualisation
......@@ -202,17 +229,18 @@ while do.epoch == 1;
end
% we may want to adopt an automatized algorithm instead of this thing
EEG = pop_eegfiltnew(EEG, 0.1, 40, 16500, 0, [], 0); % filtering 1htz lower cutof and 90upper cutof (for erp 30-40)
EEG = pop_eegfiltnew(EEG, 1, 40, 1650, 0, [], 0); % filtering 1htz lower cutof and 90upper cutof (for erp 30-40)
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
EEG = pop_epoch(EEG, {'PIC1' 'PIC3'}, [-1 3], 'epochinfo', 'yes'); % extract epochs 1 second before the picture and 3 after the picture. We may want to shortening
EEG = pop_epoch(EEG, {'PIC1' 'PIC3'}, [-0.2 1.5], 'epochinfo', 'yes'); % extract epochs 1 second before the picture and 3 after the picture. We may want to shortening
EEG = pop_rmbase(EEG, [-200 0]); % remove baseline - we establish 200 ms before the picture as the baselie.
[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET, 'overwrite', 'on'); % save
EEG = pop_subcomp(EEG, find(EEG.reject.gcompreject), 0); % remove the IC components marked as artifacts - we are statistically smoothing whatever the ICA found.
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
EEG = rejSegChan(EEG, {'EOG1','EOG2','EOG3','EOG4'}, 100); % call an algorithm for rejecting artifactual epochs and channels (see help rejSegChan) - allows a deviation of 100 micro volts before rejecting
% detect artifacts while ignoring a bunch of channels
EEG = rejSegChan(EEG, {'VEOG1','VEOG2','HEOG1','HEOG2','O9', 'O10', 'F9', 'F10', 'Fn'}, 100); % call an algorithm for rejecting artifactual epochs and channels (see help rejSegChan) - allows a deviation of 100 micro volts before rejecting
%however, if there is one channel that doens't work we don't want
%to remove all the data. We therefore check if this is a crucial
%channel and keep the data.
......@@ -230,11 +258,11 @@ while do.build == 1;
epochlist = dir([epochdir,'*.set']);
% RESEARCHER INPUT
d.photoType = {'1' '2'}; % experimental design structure - the purpose is to have a representation of the design which will be populated. This should be a minimal cutof
d.phase = {'1' '2'};
d.condition = {'1' '2' '3' '4'};
minEpoch = 12; % minimum number of epochs that needs to be retained per design ceafter artifact rejection before the subject is excluded
roi = {{'Pz' 'Cz' 'P3' 'P4'}}; % define region(s) of interest, if applicable
d.photoType = {'emo' 'ntr'};
d.phase = {'indiv' 'social'};
d.condition = {'low' 'high' 'same' 'none'};
minEpoch = 8; % minimum number of epochs that needs to be retained per design ceafter artifact rejection before the subject is excluded
roi = {}; % {'Pz' 'Cz' 'P3' 'P4'} define region(s) of interest, if applicable
% INITIALIZE VARIABLES
p = struct('times', {[]}, 'chanlocs', readlocs([scriptDir,'locfile.ced']), 'data', {{}}, 'subs', {{}}); % p is a custom structure we'll build to be used later as input to visualizaton functions
......@@ -266,7 +294,13 @@ while do.build == 1;
bkg = getBkg(EEG, d, bkg, minEpoch, readlocs([scriptDir,'locfile.ced'])); % call a function to build meta-data matrix
if bkg{end,end} == 0 % if any of the design cells has less than required trials...
s = s - 1; % ... drop the subject from the p structure
continue
end
epoch = epochList(EEG); % get information about each epoch in correct format
p.epochFields = fields(epoch)';
EEG = pop_interp(EEG, readlocs([scriptDir,'locfile.ced']), 'spherical'); % interpolate missing channels - if you threw out a channel the system will try to fill in the data
[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);
......@@ -286,15 +320,14 @@ while do.build == 1;
for condition = 1:4
t.c = strcmp(d.condition{condition},{epoch.condition});
inds = find(t.a.*t.b.*t.c); % find trials belonging to one design cell...
n(end+1) = length(inds); % ... and count them
p.data{photoType,phase,condition}(:,:,s) = mean(EEG.data(:,:,inds),3); % average trials in the selected cell
p.raw.dat{photoType,phase,condition,s} = EEG.data(:,:,inds); % trial data
p.raw.inf{photoType,phase,condition,s} = squeeze(struct2cell(epoch(inds)))';
end
end
end
if min(n) < minEpoch; % if any of the cells has less than required trials...
s = s - 1; % ... drop the subject from the p structure
end
end
save([epochdir date '-data.mat'], 'p'); % save the p structure to disk
xlswrite([epochdir date '-export.xlsx'], bkg, 'bkg'); % save the meta-data file
......@@ -303,29 +336,35 @@ end
%% VISUALIZE AND EXPORT %%
while do.analyze == 1;
buildlist = dir([epochdir,'*.mat']); % find .mat files from epochdir
load([epochdir buildlist(end).name]); % open the last one
% there are 120 emo and 60 ntr pic, therefore the maximus should be 30
% and 15, so something is off in the
% VIUSAL STATISTICS %
% plotERPeffect(data structure, factor-to-plot, cells-to-include, p-level, correctoon, time-to-plot, baseline, plot type)
% see help plotERPeffect
% Envelope plot
plotERPeffect(p,1,{1:2,1:2,1:4}, 1, 'none', [-1000 3000], [-200 0], 'envp'); % plot grand average waveforms from each channel on a single plot
plotERPeffect(p,1,{1:2,1:2,1:4}, 1, 'none', [-200 1500], [-200 0], 'envp'); % plot grand average waveforms from each channel on a single plot
% ERP plots
plotERPeffect(p,1,{1:2,1:2,1:4}, 0.05, 'fdr', [-200 1500], [-200 0], 'erp'); % e.g. photoType main effect with False Discovery rate correction
plotERPeffect(p,3,{1:2,2,1:4}, 0.05, 'cluster', [-200 1500], [-200 0], 'erp'); % e.g. condition main effect from the second phase with cluster-based statistics
% ERP plots
plotERPeffect(p,3,{1,1,[1 2]}, 0.05, 'none', [-200 1500], [-200 0], 'erp'); % e.g. photoType main effect with False Discovery rate correction
plotERPeffect(p,[3 2],{1,1:2,[1:3]}, 0.05, 'none', [-200 1500], [-200 0], 'erp'); % e.g. photoType main effect with False Discovery rate correction
plotERPeffect(p,1,{1:2,2,1:3}, 0.05, 'cluster', [-200 1500], [-200 0], 'erp'); % e.g. condition main effect from the second phase with cluster-based statistics
% SCALP PLOTS
plotERPeffect(p,1,{1:2,1:2,1:3}, 0.05, 'no', [200 300], [-200 0], 'scalpDiff'); % plot the scalp distribution of the difference between the last and first level of the selected dimension
plotERPeffect(p,3,{1:2,1:2,1:3}, 1.0, 'no', [300 400], [-200 0], 'scalpComp'); % plot the scalp distribution of the grand average
plotERPeffect(p,1,{1:2,2,[1:3]}, 0.99, 'no', [270 470], [-200 0], 'scalpDiff'); % plot the scalp distribution of the difference between the last and first level of the selected dimension
plotERPeffect(p,1,{1:2,1,1:3}, 1.0, 'no', [170 260], [-200 0], 'scalpComp'); % plot the scalp distribution of the grand average
% EXPORT (see help exportERPeffect)
exportERPeffect([epochdir date '-theta.xlsx'], p,{1:2,1:2,1:4},...
{'LPP', {'AF3' 'AF4' 'Fz' 'F3' 'F4'}, [300 800]}, [-200 0]);
exportERPeffect([epochdir date '-export.xlsx'], p,{1:2,1:2,1:4},...
{'P3', {'Pz'}, [270 470];...
'LPP', {'Pz'}, [470 800];...
'SLW' , {'Pz'}, [800 1100]}, [-200 0], 'ch', NaN, 'no' ,'yes');
do.analyze = 0;
end
FidNz 0 9.071585155 -2.359754454
FidT9 -6.711765 0.040402876 -3.251600355
FidT10 6.711765 0.040402876 -3.251600355
Fp1 -2.695405558 8.884820317 1.088308144
Fp2 2.695405558 8.884820317 1.088308144
F3 -4.459387187 6.021159964 4.365321482
F4 4.459387187 6.021159964 4.365321482
C3 -5.47913021 0.284948655 6.38332782
C4 5.47913021 0.284948655 6.38332782
P3 -5.831241498 -4.494821698 4.955347697
P4 5.831241498 -4.494821698 4.955347697
O1 -2.738838019 -8.607966849 0.239368223
O2 2.738838019 -8.607966849 0.239368223
F7 -6.399087198 4.127248875 -0.356852241
F8 6.399087198 4.127248875 -0.356852241
T3 -7.304625099 -1.866238006 -0.629182006
T4 7.304625099 -1.866238006 -0.629182006
T5 -6.034746843 -5.755782196 0.051843011
T6 6.034746843 -5.755782196 0.051843011
Fz 0 7.96264703 5.044718001
Cz 0 9.271139705 -2.211516434
Pz 0 -6.676694032 6.465208258
Oz 0 -8.996686498 0.487952047
E21 -6.518995129 2.417299399 -5.253637073
E22 6.518995129 2.417299399 -5.253637073
LM -6.174969392 -2.458138877 -5.637380998
RM 6.174969392 -2.458138877 -5.637380998
E25 -3.784983913 -6.401014415 -5.260040689
E26 3.784983913 -6.401014415 -5.260040689
E27 0 9.087440894 1.333345013
FCz 0 3.806770224 7.891304964
LE -3.743504949 6.649204911 -6.530243068
RE 3.743504949 6.649204911 -6.530243068
E31 -6.118458137 4.523870113 -4.409174427
E32 6.118458137 4.523870113 -4.409174427
Cz 0 0 8.899186843
Fp1 E1 -2.695405558 8.884820317 1.088308144
Fp2 E2 2.695405558 8.884820317 1.088308144
F3 E3 -4.459387187 6.021159964 4.365321482
F4 E4 4.459387187 6.021159964 4.365321482
C3 E5 -5.47913021 0.284948655 6.38332782
C4 E6 5.47913021 0.284948655 6.38332782
P3 E7 -5.831241498 -4.494821698 4.955347697
P4 E8 5.831241498 -4.494821698 4.955347697
O1 E9 -2.738838019 -8.607966849 0.239368223
O2 E10 2.738838019 -8.607966849 0.239368223
F7 E11 -6.399087198 4.127248875 -0.356852241
F8 E12 6.399087198 4.127248875 -0.356852241
T3 E13 -7.304625099 -1.866238006 -0.629182006
T4 E14 7.304625099 -1.866238006 -0.629182006
T5 E15 -6.034746843 -5.755782196 0.051843011
T6 E16 6.034746843 -5.755782196 0.051843011
Fz E17 0 7.96264703 5.044718001
Cz E18 0 0 8.899186843
Pz E19 0 -6.676694032 6.465208258
Oz E20 0 -8.996686498 0.487952047
F9 E21 -6.518995129 2.417299399 -5.253637073
F10 E22 6.518995129 2.417299399 -5.253637073
LM E23 -6.174969392 -2.458138877 -5.637380998
RM E24 6.174969392 -2.458138877 -5.637380998
O9 E25 -3.784983913 -6.401014415 -5.260040689
O10 E26 3.784983913 -6.401014415 -5.260040689
Fpz E27 0 9.087440894 1.333345013
FCz E28 0 3.806770224 7.891304964
VEOG1 LE -3.743504949 6.649204911 -6.530243068
VEOG2 RE 3.743504949 6.649204911 -6.530243068
HEOG1 E31 -6.118458137 4.523870113 -4.409174427
HEOG2 E32 6.118458137 4.523870113 -4.409174427
......@@ -23,17 +23,20 @@ pl1way <- function(dat,y,x) {
```
###some dataset definision
###social - second phase
###individual - first phase.
```{r get and convert the file, include=FALSE}
setwd("C:/Users/Amit/Dropbox/Research/emotional conformity-non conformity/motivated_contation_EEG_2017/Results")
#setwd("C:/Users/ro/OneDrive - Tartu likool/Teadus/REGULATION/16 motivated contagion")
d <- read.csv("eeg_data_motivated_contagion.csv")
d=d[,c(1,3:5,9,11,14,17:18,21:22,24:26)]
d$grate = ifelse(d$grate == 0, NA, d$grate)
head(d)
#d = subset (d, condition != "none")
#ds = subset (d, phase == "social")
##remove participants
d = subset (d, subject != "5787127" & subject != "5858527")
# separate phases
a <- d[d$phase=="indiv",]
......@@ -47,35 +50,66 @@ colnames(d1) <- gsub(".x", ".indiv", colnames(d1), fixed =T)
colnames(d1) <- gsub(".y", ".social", colnames(d1), fixed =T)
# remove "none"" and fix comparison levels
d1 <- d1[d1$condition!="none",]
head(d1)
d1e = subset (d1,photoType =="emo" )
d1n = subset (d1,photoType =="ntr" )
#d1 <- d1[d1$condition!="none",]
#I just want to make sure that this just organizes the condition
d1$condition <- factor(as.character(d1$condition), levels = c( "low", "high","same"))
d1$condition <- factor(as.character(d1$condition), levels = c( "low", "high","same", "none"))
d1$photoType <- relevel(d1$photoType,"ntr")
# add phase differences
d1$P3.dif = d1$P3.indiv -d1$P3.social
d1$LPP.dif = d1$LPP.indiv- d1$LPP.social
d1$SLW.dif = d1$SLW.indiv-d1$SLW.social
d1$rate.dif = d1$rate.indiv -d1$rate.social
d1$P3.dif = d1$P3.social-d1$P3.indiv; hist (d1$P3.dif)
d1$LPP.dif = d1$LPP.social-d1$LPP.indiv
d1$SLW.dif = d1$SLW.social-d1$SLW.indiv
d1$rate.dif = d1$rate.social-d1$rate.indiv
# add difference between group and individual ratings
d1$rateDifCon <- d1$grate - d1$rate.social
d1$rateDifCon <- d1$rate.social-d1$grate
d1e = subset (d1,photoType =="emo" &(rate.indiv<8 &rate.indiv ))
d1n = subset (d1,photoType =="ntr" )
d1$photoType <- relevel(d1$photoType,"ntr")
```
### perliminary analysis, get raw means
```{r}
means = d1e%>%
group_by(photoType,condition) %>%
summarise(
P3s =mean(P3.social, na.rm = T),
LPPs = mean(LPP.social, na.rm = T),
SLWs = mean(SLW.social , na.rm =T),
P3i = mean(P3.indiv, na.rm = T),
LPPi = mean(LPP.indiv, na.rm = T),
SLWi = mean(SLW.indiv, na.rm = T),
P3.dif = mean(P3.dif, na.rm = T),
LPP.dif = mean(LPP.dif, na.rm = T),
SLW.dif = mean(SLW.dif, na.rm = T)
)
data.frame(means)
rm (means)
```
### A couple of descriptive pictures
```{r descriptives, warning=F}
# frequency table of trials
table(d1$subject, d1$condition, d1$photoType)
a=table(d1$subject, d1$condition)
length(table(mean(data.frame(apply(a,1,sum))[,1])))
table(d1e$subject, d1e$condition)
# independent variables
ggpairs(d1e[,c("rate.social","grate", "rateDifCon")])
......@@ -88,10 +122,10 @@ ggpairs(d1[,c("P3.indiv","P3.social", "P3.dif")])
ggpairs(d1[,c("LPP.indiv","LPP.social", "LPP.dif")])
ggpairs(d1[,c("SLW.indiv","SLW.social", "SLW.dif")])
ggpairs(d1[,c("P3.indiv","P3.social")])
ggpairs(d1[,c("LPP.indiv","LPP.social")])
ggpairs(d1e[,c("P3.indiv","P3.social")])
ggpairs(d1e[,c("rate.indiv","LPP.social")])
ggpairs(d1[,c("SLW.indiv","SLW.social")])
ggpairs(d1[,c("rate.indiv","rate.social")])
ggpairs(d1e[,c("rate.indiv","rate.social")])
# Do we need to control for order? No
cor(d[,c("P3","LPP", "SLW", "rate", "grate")],d$order, use="pairwise.complete.obs")
......@@ -111,13 +145,21 @@ The brain-centric model assumes that the first presentation response (P3.indiv t
The alternative is of course to use what we manipulated to regress out brain variance not related to the social manipulation. If we do that, we shouldn't any more include the first phase EEG, because we then run the risk of regressing the same thing out twice.
### 2) For brain-centric approach, should we use a difference score or residualization?
In principle using the difference between P3.social and P3.indiv as a depent variable is also a way to do the brain-centric approach. There is some reason to think residualization works better, because variance from the first phase will "leak" into difference scores but not into residualized scores:
In principle using the difference between P3.social and P.indiv as a depent variable is also a way to do the brain-centric approach. There is some reason to think residualization works better, because variance from the first phase will "leak" into difference scores but not into residualized scores:
```{r individual and social , warning=F}
head(d1)
rsd <- lmer(P3.social ~ P3.indiv + (1|subject), data = d1) ; summary (rsd)
r = lmer (P3.social ~ rate.indiv +(1|subject) , d1e); summary(r)
r = lmer (P3.social ~ rate.social+grate+(1|subject) , d1e); summary(r)
r = lmer (LPP.social ~ grate+rate.social +(1|subject) +(1|photo), d1e); summary(r)
r = lmer (SLW.social ~ rate.indiv +(1|subject) , d1e); summary(r)
r = lmer (SLW.social ~ rate.social+grate +(1|subject) , d1e); summary(r)
#rsd <- lmer(P3.social ~ P3.indiv + (1|subject), data = d1e) ; summary (rsd)
ggpairs(cbind(d1[,c("P3.indiv", "P3.social","P3.dif")],P3.residuals=resid(rsd)))
......@@ -136,11 +178,14 @@ One reaosn to prefer the factor is that this allows us to detect nonlinear effec
Based on various reasons, including some saturday-night-p-hacking, I stuck to the brain-centric models using group ratings as factor:
### P3
### P3 analysis
```{r P3 dif, warning=F}
# simplest model
summary(fit <- lmer(P3.social ~ condition +P3.indiv+order.indiv+order.social + (1|subject), data = d1))
summary(fit <- lmer(P3.social ~ condition +grate+P3.indiv+order.indiv+order.social+ (1|subject)+(1|photo), data = d1e))
# get p-values for the random effects
rand(fit)
......@@ -162,7 +207,27 @@ pl1way(cbind(D3,dep=resid(rsd)),"dep", "condition")
```{r, warning=F}
# simplest model
summary(fit <- lmer(LPP.social ~ LPP.indiv + condition + (1|subject), data = d1))
summary(fit <- lmer(LPP.social ~ condition +LPP.indiv+ (1|subject), data = d1e))
270-470
470-800
200/530
330/530
d1e$a = d1e$LPP.social*0.3773585
d1e$b = d1e$LPP.social*0.6226415
d1e$c = d1e$LPP.indiv*0.3773585
d1e$d = d1e$LPP.indiv*0.6226415
d1e$early.social= (rowMeans(d1e[,c('a', 'b')], na.rm=T))
d1e$early.indiv= (rowMeans(d1e[,c('c', 'd')], na.rm=T))
summary(fit <- lmer(early.social ~condition +early.indiv+(1|photo)+(1|subject), data = d1e))
# get p-values for the random effects
rand(fit)
......@@ -184,9 +249,9 @@ pl1way(cbind(D3,dep=resid(rsd)),"dep", "condition")
```{r, warning=F}
# simplest model
summary(fit <- lmer(SLW.social ~ SLW.indiv + condition + (1|subject), data = d1))
summary(fit <- lmer(SLW.social ~ condition +SLW.indiv+ order.indiv+order.social+ (1|subject)+(1|photo), data = d1e))
# get p-values for the random effects
# get p-values for the random effects
rand(fit)
# get a significance test for all terms
anova(fit)
......@@ -209,23 +274,92 @@ Based on these analyses, we have a small but significant effect on P3 and no eff
```{r dif, include=FALSE}
summary(fit <- lmer(scale(P3.dif, scale = F) ~ condition+ (1|subject) + (1|photo), data = d1[d1$photoType=="emo",]))
summary(fit <- lmer(P3.social ~ condition +P3.indiv+order.indiv+order.social + (1|subject), data = d1))
summary(fit <- lmer(scale(P3.dif, scale = F) ~ condition+ (1|subject) + (1|photo), data = d1e))
summary(fit <- lmer(P3.indiv ~ condition + (1|subject), data = d1))
summary(fit <- lmer(scale(LPP.dif, scale = F) ~ condition +order.social+order.indiv+ (1|subject) + (1|photo), d1e))
summary(fit <- lmer(scale(SLW.dif, scale = F) ~ condition + (1|subject) + (1|photo), data = d1[d1$photoType=="emo",]))
```
summary(fit <- lmer(scale(LPP.dif, scale = F) ~ condition + (1|subject) + (1|photo), data = d1[d1$photoType=="emo",]))
```{r order}
summary(fit <- lmer(LPP.social ~ order.social+ (1|subject) + (1|photo)-1, d1e))
summary(fit <- lmer(LPP.indiv ~ order.indiv+ (1|subject) + (1|photo)-1, d1e))
```
###actual ratings
```{r actual ratings}
r = lmer(rate.indiv~ condition+(1|subject) +(1|photo),d1e);summary (r)
r = lmer(scale(rate.dif, scale =F)~ condition+order.indiv+order.social -1+(1|subject) +(1|photo),d1e);summary (r)
head(d1e)
r = lmer (rate.social ~ condition +order.indiv+order.social + (1|subject) + (1|photo), d1e) ; summary (r)
```
####reaction time
```{r}
r = lmer(rt.social~ condition+order.indiv+order.social +(1|subject) +(1|photo),d1e);summary (r)
##correlate reaction time with lpp
r = lmer(LPP.social~ rt.social+order.indiv+order.social +(1|subject) +(1|photo),d1e);summary (r)
```
### analysis of individual participants
```{r}
addmargins(table(d1e$subject, d1e$condition))
models = d1e%>%
group_by(subject) %>%
do(r = print(summary(lm(LPP.social ~ condition +P3.indiv, data = .))))
coef = models %>%
do(data.frame (low = (.$r)$coefficients[[1]],
high = (.$r)$coefficients[[2]],
same = (.$r)$coefficients[[1]]
))
coef =cbind(models[,1], coef)
coef$test= ifelse(coef$high>1,1,2)
table(coef$test)
coef=arrange (coef,coef[,1])
rm(models)
summary(fit <- lmer(scale(SLW.dif, scale = F) ~ condition + (1|subject) + (1|photo), data = d1[d1$photoType=="emo",]))
```
###compute differences from neutral
```{r compute difference from neutral}
head(d)
##creeat dataset of the averages.
......@@ -269,46 +403,4 @@ r= lm (SLW~ condition+SLWi.e+SLWi.n, dd); summary (r)
```
###actual ratings
```{r actual ratings}
d2 = subset (d1, rate.social <8 &rate.social>1)
r = lmer(scale(rate.dif, scale =F)~ condition+order.indiv -1+(1|subject) +(1|photo),d2);summary (r)
```
```{r analysis slw, include=FALSE}
head(pl1way)
# AFFECTIVE MAIN EFFECT
summary(fit <- lmer(SLW.indiv ~ photoType + (1|subject) + (1|photo), data = d1))
pl1way(D3,"P3.indiv","photoType")
pl1way(D3,"LPP.indiv","photoType")
pl1way(D3,"LPP.indiv","photoType")
# continuous model
summary(fit <- lmer(SLW.social ~ SLW.indiv + rate*grate + (1|subject) + (1|photo), data = D3))
coefplot(fit)
visreg(fit, xvar = "grate", by = "rate")
# categorical model
summary(fit <- lmer(SLW.social ~ SLW.indiv + photoType*condition + (1|subject) + (1|photo), data = D3))
coefplot(fit)
# model without neutral pictures
summary(fit <- lmer(SLW.social ~ SLW.indiv + condition + (1|subject) + (1|photo), data = D3[D3$photoType=="emo",]))
coefplot(fit)
```
function epoch = epochList(EEG)
epoch = struct(); % genereerib trialitestruktuuri, kus oleks vaid ks vrtus triali kohta
epoch = struct(); % genereerib trialitestruktuuri, kus oleks vaid ks vrtus triali kohta
for e = 1:EEG.trials
fld = fields(EEG.epoch)';
for f = 2:length(fld)
if iscell(EEG.epoch(e).(fld{f})) % kui on hes trialis mitu eventi
EEG.epoch(e).(fld{f}) = EEG.epoch(e).(fld{f}){end};
if iscell(EEG.epoch(e).(fld{f})) % kui on hes trialis mitu eventi
EEG.epoch(e).(fld{f}) = EEG.epoch(e).(fld{f}){1};
end
if ischar(EEG.epoch(e).(fld{f}))
epoch(e).(fld{f}(6:end)) = strtrim(EEG.epoch(e).(fld{f})); % kui on string, eemalda thikud
epoch(e).(fld{f}(6:end)) = strtrim(EEG.epoch(e).(fld{f})); % kui on string, eemalda thikud
else
% epoch(e).(fld{f}(6:end)) = EEG.epoch(e).(fld{f});
epoch(e).(fld{f}(6:end)) = num2str(EEG.epoch(e).(fld{f}));
......
......@@ -121,7 +121,7 @@ function outf = exportERPeffect(savefile, inp, design, comp, baseline, loc, late
for d = 1:numel(data) % loop through experimental conditions