############################################### ### Mass univariate statistics for ERP data ### ############################################### # Written by Emily Coderre, PhD, 2015-10-26 # This function runs mass univariate statistics on a given array of EEG data in R. # The options for methods are based on those described in Groppe, Urbach, & Kutas (2011): Mass univariate analysis of event-related brain potentials/fields I: A critical tutorial review. Psychophysiology, 48, 1711-1725. Refer to this paper for an overview of the pros and cons of each method and to determine which method is best to use in different situations. # The function takes an input array of EEG data of the format [subjects conditions timepoints electrodes] and outputs a matrix of [electrodes timepoints] with a significance mask of 0 (not significant) or 1 (significant) according to whether every point estimate is significant based on the mass univariate method of choice. This significance mask can then be used to plot the results of the mass univariate tests massunistats <- function(alldata,method,submethod,electrodes = c(1:dim(alldata)[4]),timepoints = c(1:dim(alldata)[3]), alpha = 0.05, nIterations = 5000, successive = 2, u = 1,tails=2){ ### Arguments to function massunistats: # - alldata: four-dimensional matrix of [subjects conditions timepoints electrodes]. **Currently only works with two conditions and with within-subject observations # - method: # 'perm' = permutation # 'FDR' = false discovery rate # 'GFWER' = generalized family-wise error rate # - submethod: if method is 'perm' or 'FDR', specify submethod # - method 'perm': # submethod 'tmax' = strong control of FWER via t-max statistic # submethod 'cluster' = cluster-based permutation with weak control of FWER. Currently uses cluster mass statistic. # - method 'FDR' # submethod 'BH' = Benjamini & Hochberg method # submethod 'BY' = Benjamini & Yekutieli method # submethod 'BKY' = Benjamini, Krieger & Yekutieli method # - If method 'GFWER', no submethod is needed # - electrodes (optional): specify which electrodes to look at for the analyses. Input should be numerical, corresponding to the electrode number (e.g. Pz = electrode #101 on an EGI Hydrocel 256-channel system). Default: use all electrodes (takes the number of electrodes from the 4th dimension of 'alldata'). Groppe et al. (2011) recommend limiting analyses to electrodes that are expected to show effects, to reduce the number of tests # - timepoints (optional): specify which timepoints to look at for the analyses. Input should correspond to the recorded timepoint, not the time in ms. For example, 100 ms after stimulus onset would correspond to timepoint #26 at a sampling rate of 250Hz (e.g. seq(0,100,4)). Default: use all timepoints (takes the number of timepoints from the 3rd dimension of 'alldata'). Groppe et al (2011) suggest downsampling to about 100 Hz to reduce the number of tests. # - alpha (optional): alpha level for significance tests. Default: 0.05. # - nIterations (optional): number of iterations to run for permutation tests. Default is 5000. Groeppe et al. (2011) and Manly (1997) recommend 1000 permutations for 5% alpha level, 5000 for 1% alpha. This can estimate p-values within 2%. More (e.g. 10,000) can be run if better precision is required. # - successive (optional): number of successive time windows required to define a cluster (for cluster-based permutation only). Default: 2 (based on the parameter used in an example given in Groppe et al. 2011). # - u (optional): threshold for number of false discoveries (for GFWER only). By allowing a small number of false discoveries, the power of the analysis increases over FWER. If u = 0, then GFWER is equal to strong control of FWER. Default: (based on the parameter used in an example given in Groppe et al. 2011). # - tails (optional): one-tailed or two-tailed tests. Default: 2. ### Output: # - sigmask: array of [electrodes timepoints] indicating whether the test was significant (1) or not significant (0) at each point estimate. # Require submethod argument for 'perm' and 'FDR' if ((method == 'perm' | method == 'FDR') & missing(submethod)){ stop("Please specify submethod") } # Summarize the method and print to console cat(paste('Running mass univariate stats using method "',method,'", submethod "',submethod,'"\n',sep='')) ## Select electrodes and timepoints from alldata ## data = alldata[,,timepoints,electrodes] ################################# ### Permutation-based methods ### ################################# if (method == 'perm'){ ### Calculate observed test statistic for every timepoint/electrode of interest obsTs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ CONmatrix = data[,1,,e] INCONmatrix = data[,2,,e] for (t in 1:length(timepoints)){ obsTs[e,t] = t.test(CONmatrix[,t],INCONmatrix[,t],paired=TRUE)$stat } } ############################################ ### Strong control of FWER: t-max method ### if (submethod == 'tmax'){ ## Run permutations ## nullDist = vector() for (it in 1:nIterations){ cat(paste('Thinking: Iteration ',it,'\n',sep='')) ### Permute conditions within subject for every timepoint and electrode ### # Recommended to randomize conditions within subject if comparison is within-subject # Recommended to randomize subject between conditions if comparison is between-subject randdata = array(dim=dim(data)) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ for (s in 1:dim(data)[1]){ # subject loop subdata = data[s,,t,e] # take that subject's data for each condition at the specified timepoint and electrode randdata[s,,t,e] = sample(subdata) # randomize the conditions within each subject } } } ### Get matrix of permuted test statistics for each timepoint and electrode in permuted data (randdata) ### simTs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ xCONmatrix = randdata[,1,,e] xINCONmatrix = randdata[,2,,e] for (t in 1:length(timepoints)){ simTs[e,t] = t.test(xCONmatrix[,t],xINCONmatrix[,t],paired=TRUE)$stat } } # Find most extreme test statistic over entire family of tests (i.e. all electrodes/timepoints) nullDist[it] = max(abs(simTs)) } # iterations ### Derive p-value from null distribution of permuted scores if (tails == 1){quantileCutoff = 1-alpha} # one-tailed if (tails == 2){quantileCutoff = 1-(alpha/2)} # two-tailed # Get the critical t-value for the simulated distribution critT = quantile(nullDist,quantileCutoff) ### For each timepoint/electrode, determine whether observed t-value is greater than critical t-value # Create a significance mask of 0s (ns) and 1s (sig) for every timepoint and electrode # sigmask is [electrodes timepoints] sigmask = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ # time loop if (obsTs[e,t] > critT){sigmask[e,t] = 1} else {sigmask[e,t] = 0} } } } # tmax method ########################################## ### Cluster-based weak control of FWER ### if (submethod == 'cluster'){ ## Run permutations ## nullDist = vector() for (it in 1:nIterations){ cat(paste('Thinking: Iteration ',it,'\n',sep='')) ### Permute conditions within subject for every timepoint and electrode ### # Recommended to randomize conditions within subject if comparison is within-subject # Recommended to randomize subject between conditions if comparison is between-subject randdata = array(dim=dim(data)) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ for (s in 1:dim(data)[1]){ # subject loop subdata = data[s,,t,e] # take that subject's data for each condition at the specified timepoint and electrode randdata[s,,t,e] = sample(subdata) # randomize the conditions within each subject } } } ### Get matrix of permuted test statistics for each timepoint and electrode in permuted data (randdata) ### simTs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ xCONmatrix = randdata[,1,,e] xINCONmatrix = randdata[,2,,e] for (t in 1:length(timepoints)){ simTs[e,t] = t.test(xCONmatrix[,t],xINCONmatrix[,t],paired=TRUE)$stat } } # 1) From uncorrected t-distribution (simTs), define a threshold (e.g. p<0.05) and throw out any t-values that do not exceed the threshold # Find critT corresponding to p<0.05 from uncorrected distribution if (tails == 1){quantileCutoff = 1-alpha} # one-tailed if (tails == 2){quantileCutoff = 1-(alpha/2)} # two-tailed critTuncorr = quantile(simTs,quantileCutoff) # 2) Group significant t-scores into clusters at adjacent timepoints/electrodes found = 0 sig = 0 notsig = 0 onset = vector() offset = vector() cnt = 1 successive = 2 clusterTs = vector() clusterTcount = 1 for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ ## Determine whether t-value is greater than threshold. Those below the uncorrected critT threshold are not considered as being in clusters if (abs(simTs[e,t]) > critTuncorr){ # if t-value is above threshold, add 1 to sig sig = sig + 1 if (sig >= successive & found == 0){ # if number of sig windows has reached successive threshold, mark that as an onset onset[cnt] = t - (successive-1) found = 1 # mark as found } } if (abs(simTs[e,t]) < critTuncorr & found == 0){ # if t-value is not above threshold and was not before, return sig to 0 sig = 0 } if (abs(simTs[e,t]) < critTuncorr & found == 1){ # if t-value is not above threshold but previously was, return sig to zero and mark offset sig = 0 found = 0 offset[cnt] = t-1 cnt = cnt+1 } } # timepoints # if significance continues to the end of the window, offset will be shorter than onset, so just fill in the last simTspoint if (length(onset) == length(offset)+1){ offset[cnt] = dim(simTs)[2] } ### 3) t-scores of each cluster are summed to produce a cluster-level t-score, i.e. "mass" of the cluster # creates a vector of all clustered t-score masses over all electrodes if (length(onset)>0){ for (cls in 1:length(onset)){ clusterTs[clusterTcount] = sum(abs(simTs[e,onset[cls]:offset[cls]])) clusterTcount = clusterTcount + 1 } } ### 4) Most extreme cluster-level t-score across all timepoints and electrodes is used to derive a null distribution if (length(clusterTs) > 0){ nullDist[it] = max(clusterTs) } else {nullDist[it] = 0} } # electrodes } # iterations ### Derive critical cluster mass value from null distribution of permuted scores if (tails == 1){quantileCutoff = 1-alpha} # one-tailed if (tails == 2){quantileCutoff = 1-(alpha/2)} # two-tailed # Get the critical test statistic for the simulated distribution critT = quantile(nullDist,quantileCutoff) ### Calculate the mass of each cluster in the observed data and compare to the null distribution of cluster masses # Find critT corresponding to p<0.05 from uncorrected distribution if (tails == 1){quantileCutoff = 1-alpha} # one-tailed if (tails == 2){quantileCutoff = 1-(alpha/2)} # two-tailed critTuncorr = quantile(obsTs,quantileCutoff) # Group significant t-scores into clusters at adjacent timepoints/electrodes # Note: this currently only goups by temporal clusters, not spatial clusters. Groppe et al (2011) define an electrode's "spatial neighborhood" as all electrodes within 5.4 cm. If looking at only a select few electrodes distributed across the scalp, spatial clustering is likely not needed; however, for exploratory analyses looking at many adjacent electrodes, you will need to implement a spatial cluster in the code below. # Create a significance mask of 0s (ns) and 1s (sig) for every timepoint and electrode # sigmask is [electrodes timepoints] sigmask = matrix(0,length(electrodes),length(timepoints)) # make sigmask a matrix of zeros - this will essentially make everything ns by default, and then add in significant clusters found = 0 sig = 0 notsig = 0 onset = vector() offset = vector() cnt = 1 successive = 2 obsclusterTs = vector() clusterTcount = 1 for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ ## Determine whether t-value is greater than threshold. Those below the uncorrected critT threshold are not considered as being in clusters if (abs(obsTs[e,t]) > critTuncorr){ # if t-value is above threshold, add 1 to sig sig = sig + 1 if (sig >= successive & found == 0){ # if number of sig windows has reached successive threshold, mark that as an onset onset[cnt] = t - (successive-1) found = 1 # mark as found } } if (abs(obsTs[e,t]) < critTuncorr & found == 0){ # if t-value is not above threshold and was not before, return sig to 0 sig = 0 } if (abs(obsTs[e,t]) < critTuncorr & found == 1){ # if t-value is not above threshold but previously was, return sig to zero and mark offset sig = 0 found = 0 offset[cnt] = t-1 cnt = cnt+1 } } # timepoints # if significance continues to the end of the window, offset will be shorter than onset, so just fill in the last simTspoint if (length(onset) == length(offset)+1){ offset[cnt] = length(obsTs) } ### For each cluster, determine whether the mass of that cluster is larger than the critical test statistic. if (length(onset)>0){ for (cls in 1:length(onset)){ obsclusterT = sum(abs(obsTs[e,onset[cls]:offset[cls]])) # get mass of cluster # if observed cluster mass exceeds critical cluster mass, mark as significant if (obsclusterT > critT){ sigmask[e,onset[cls]:offset[cls]] = 1 } else {sigmask[e,onset[cls]:offset[cls]] = 0} } } else {sigmask[e,] = 0} # if no clusters found in observed data, mark everything as ns } # electrodes } # cluster submethod } # perm method #################################### ### False Discovery Rate control ### #################################### if (method == 'FDR'){ ## Select which method to use # Method 1 = Benjamini & Hochberg (BH) # Method 2 = Benamini & Yekutieli (BY) # Method 3 = Benjamini, Krieger, and Yekutieli (BKY 2006) ##################################### # Calculate p-values at every timepoint/electrode obsPs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ CONmatrix = data[,1,,e] INCONmatrix = data[,2,,e] for (t in 1:length(timepoints)){ obsPs[e,t] = t.test(CONmatrix[,t],INCONmatrix[,t],paired=TRUE)$p.value } } ######################################## # Sort p-values from entire family of m tests from smallest-largest, such that Pi is ith largest p-value allPs = sort(obsPs) m = length(allPs) # m = number of tests ########################################## ## For each i, calculate whether a specific formula is true # Create a significance mask of 0s (ns) or 1s (sig) sigmask = array(dim=c(length(electrodes),length(timepoints))) ### 3a) Benjamini and Hochberg (BH) procedure: Pi <= (i/m)*alpha if (submethod == 'BH'){ iLogical = vector() for (i in 1:m){ iLogical[i] = allPs[i] <= (i/m)*alpha } # If at least one value of i satisfies the relationship, define k as the largest value of i for which the above algorithm is true, then reject the null hypothesis for hypotheses 1:k if (length(which(iLogical == TRUE)) > 0){ # if there were any tests that were true k = max(which(iLogical == TRUE)) # Define k as the critical p-value, and everything below that will be significant critP = allPs[k] # Create significance mask of 1 (sig) or 0 (ns) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ if (obsPs[e,t] <= critP){ sigmask[e,t] = 1 } else {sigmask[e,t] = 0} } } } else {sigmask[,] = 0} # otherwise, mark all points as 0 (ns) } # BH method ### 3b) Benjamini and Yekutieli (BY) procedure: Pi <= (i/(m*(sigma1:m of 1/j)*alpha if (submethod == 'BY'){ iLogical = vector() for (i in 1:m){ # evaluate sigma by running a loop sigma = 0 for (j in 1:m){sigma = sigma + (1/j)} iLogical[i] = allPs[i] <= (i/(m * sigma)) * alpha } # If at least one value of i satisfies the relationship, define k as the largest value of i for which the above algorithm is true, then reject the null hypothesis for hypotheses 1:k sigmask = array(dim=c(length(electrodes),length(timepoints))) if (length(which(iLogical == TRUE)) > 0){ # if there were any tests that were true k = max(which(iLogical == TRUE)) # Define k as the critical p-value, and everything below that will be significant critP = allPs[k] # Create significance mask of 1 (sig) or 0 (ns) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ if (obsPs[e,t] <= critP){ sigmask[e,t] = 1 } else {sigmask[e,t] = 0} } } } else {sigmask[,] = 0} # otherwise, mark all points as 0 (ns) } # BY method ### 3c) Benjamini, Krieger, and Yekutieli (BKY 2006) procedure: Pi <= (i/m) * alphaPrime if (submethod == 'BKY'){ ## Stage 1 ## alphaPrime = alpha/(1+alpha) iLogical = vector() for (i in 1:m){ iLogical[i] = allPs[i] <= (i/m)*alphaPrime } # if there were any tests that were true if (length(which(iLogical == TRUE)) > 0){ # define k as highest value of i for which equation is true and reject H0 for tests 1:k k = max(which(iLogical == TRUE)) # define r1 as number of hypotheses rejected by Stage 1 r1 = k # if r1 == 0 or r1 == m, no hypotheses are rejected and procedure stops # else continue to Stage 2 # if all tests were true, stop procedure and mark all as significant if (k == m){ sigmask[,] = 1 } else { ## Stage 2 ## alphaDoublePrime = (m/(m-r1))*alphaPrime iLogical = vector() for (i in 1:m){ iLogical[i] = allPs[i] <= (i/m)*alphaDoublePrime } if (length(which(iLogical == TRUE)) > 0){ # if there were any tests that were true k = max(which(iLogical == TRUE)) # Define k as the critical p-value, and everything below that will be significant critP = allPs[k] # Create significance mask of 1 (sig) or 0 (ns) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ if (obsPs[e,t] <= critP){ sigmask[e,t] = 1 } else {sigmask[e,t] = 0} } } } else {sigmask[,] = 0} # otherwise, mark all points as 0 (ns) } } else {sigmask[,] = 0} # if k == 0 (no tests fulfilled equation), mark all points as 0 (ns) } # BKY method } # FDR method ################################################# ### Generalized Familywise Error Rate (GFWER) ### ################################################# # Method from Korn et al 2004 (KTMS 2004 method) # Guarantees that the number of false discoveries does not exceed a specified value u if (method == 'GFWER'){ ##################################### # Calculate p-values at every timepoint/electrode obsPs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ CONmatrix = data[,1,,e] INCONmatrix = data[,2,,e] for (t in 1:length(timepoints)){ obsPs[e,t] = t.test(CONmatrix[,t],INCONmatrix[,t],paired=TRUE)$p.value } } ######################################## # Sort p-values from entire family of m tests from smallest-largest, such that Pi is ith largest p-value allPs = sort(obsPs) m = length(allPs) # m = number of tests ##############################################3 # For all tests, obtain adjusted p-values via permutation tests. Instead of using most extreme test statistic, use (u+1)th most extreme statistic nullDist = vector() for (it in 1:nIterations){ cat(paste('Thinking: Iteration ',it,'\n',sep='')) ### Permute conditions within subject for every timepoint and electrode ### # Recommended to randomize conditions within subject if comparison is within-subject # Recommended to randomize subject between conditions if comparison is between-subject randdata = array(dim=dim(data)) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ for (s in 1:dim(data)[1]){ # subject loop subdata = data[s,,t,e] # take that subject's data for each condition at the specified timepoint and electrode randdata[s,,t,e] = sample(subdata) # randomize the conditions within each subject } } } ### Get matrix of p-values for each timepoint and electrode in permuted data (randdata) ### simPs = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ xCONmatrix = randdata[,1,,e] xINCONmatrix = randdata[,2,,e] for (t in 1:length(timepoints)){ simPs[e,t] = t.test(xCONmatrix[,t],xINCONmatrix[,t],paired=TRUE)$p.value } } ### Find (u+1)th most extreme test statistic over entire family of tests (i.e. all electrodes/timepoints) # Sort all simPs from smallest-largest and take value u+1 to put into null distribution nullDist[it] = sort(simPs)[u+1] } # iterations ### Derive p-value from null distribution of permuted scores if (tails == 1){quantileCutoff = 1-alpha} # one-tailed if (tails == 2){quantileCutoff = 1-(alpha/2)} # two-tailed # Get the critical p-value for this simulated distribution critP = quantile(nullDist,quantileCutoff) # For each timepoint/electrode, determine whether observed p-value is less than critical p-value. # Create a significance mask of 0s (ns) and 1s (sig) for every timepoint and electrode sigmask = array(dim=c(length(electrodes),length(timepoints))) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ # time loop if (obsPs[e,t] < critP){sigmask[e,t] = 1} else {sigmask[e,t] = 0} } } ################################## # Automatically reject H0 for u tests with smallest p-values (give them adjusted p-value of 0) # Note: Groppe et al say this should be done first, but it's easier to implement in programming if done after. Should give the same results critP_u = allPs[u] # Create significance mask of 1 (sig) or 0 (ns) for (e in 1:length(electrodes)){ for (t in 1:length(timepoints)){ if (obsPs[e,t] <= critP_u){ sigmask[e,t] = 1 } } } } # GFWER method return(sigmask) } ######################################### ### Plot results of permutation tests ### ######################################### ## Function to plot waveforms and results of mass univariate stats plotmassunistats <- function(alldata,electrode,electrodes = c(1:dim(alldata)[4]),timepoints = c(1:dim(alldata)[3]),condnames,plotcolors,xmin,xmax,ymin,ymax,srate){ ### Arguments to function plotmassunistats: # - alldata: four-dimensional matrix of [subjects conditions timepoints electrodes]. # - electrode: electrode number to plot # - electrodes: all electrodes included in the mass univariate statistics. Input should be numerical, corresponding to the electrode number (e.g. Pz = electrode #101 on an EGI Hydrocel 256-channel system). Default: use all electrodes (takes the number of electrodes from the 4th dimension of 'alldata'). # - timepoints (optional): the timepoints that were included in the mass univariate statistics. Input should correspond to the recorded timepoint, not the time in ms. For example, 100 ms after stimulus onset would correspond to timepoint #26 at a sampling rate of 250Hz (e.g. seq(0,100,4)). Default: use all timepoints (takes the number of timepoints from the 3rd dimension of 'alldata'). # - condnames: vector of condition names # - plotcolors: vector of colors to plot for each condition. Include 2 colors for each of your 2 conditions, plus a third color for plotting the results of the mass univariate statistics in bars along the bottom # - xmin: minimum value of x-axis to plot # - xmax: maximum value of x-axis to plot # - ymin: minimum value of y-axis to plot # - ymax: maximum value of y-axis to plot # - srate: sampling rate in Hz (e.g. 250) ### Collapse over subjects for plotting ### alltimepoints = dim(alldata)[3] allelectrodes = dim(alldata)[4] avdata = array(dim=c(length(condnames),alltimepoints,allelectrodes)) for (c in 1:length(condnames)){ for (t in 1:alltimepoints){ for (e in 1:allelectrodes){ avdata[c,t,e] = mean(alldata[,c,t,e]) } } } ### Define plotting parameters ### xa = seq(xmin,(xmax-(1000/srate)),1000/srate) # vector of timepoints ax = seq(xmin,xmax,100) # vector for tick marks ax1 = seq(labels,xmax,200) # vector for labels ay = seq(ymin,ymax,1) plot(x=c(xmin,xmax),y=c(0,0),type="l",ylim=c(ymax,ymin),xlab="",ylab="",axes=FALSE) lineWidth = 2 lines(x=c(0,0),y=c(-0.25,0.25),ylim=c(ymax,-ymin),col="black") axis(ax,side=1,labels=FALSE,pos=0) # axis for 100ms tick marks axis(ax1,side=1,labels=ax1,pos=0) # 200 ms labels axis(ay,side=2,labels=ay,pos=xmin) # y axis ### Plot waveforms ### for (c in 1:length(condnames)){ lines(x=xa,y=avdata[c,,electrode],col=plotcolors[c],lwd=3) } legend('topright',condnames, col = c(plotcolors), text.col = "black", lty = 'solid',lwd=4, bty="n",merge = TRUE) ## Draw colored bars to show results of mass univariate tests ### e = which(electrodes==electrode) for (t in 1:length(timepoints)){ if (sigmask[e,t] == 1){ # if significant segments(xa[timepoints[t]],ymax-0.5,xa[timepoints[t]],ymax-0.3,col=plotcolors[3],lwd=4) } } }